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ABSTRACT 


NASA  and  private  launch  providers  have  a  need  to  understand  and  control  ignition  overpressure 
blast  waves  that  are  generated  by  a  solid  grain  rocket  during  ignition.  Research  in  accurate  com¬ 
putational  fluid  dynamics  prediction  of  the  launch  environment  is  underway.  A  clearer  picture 
is  emerging  from  empirical  data  which  more  precisely  categorizes  all  the  dissipative  mecha¬ 
nisms  present  in  droplet-shock  interactions.  In  this  dissertation,  water  droplets  and  their  effects 
due  to  vaporization  are  represented  as  a  control  action  and  two  new  optimal  control  problems 
are  formulated  concerning  unsteady  shock  wave  attenuation.  A  single-phase  control  problem  is 
formulated  by  representing  the  effect  of  droplet  vaporization  as  an  energy  sink  on  the  right  hand 
side  of  the  unsteady  Euler  Equations  in  one  dimension.  Results  for  the  optimal  distribution  of 
equivalent  mass  of  water  vaporized  for  a  given  level  of  attenuation  are  presented.  A  two-phase 
control  problem  consists  of  solving  for  the  initial  optimal  water  droplet  distribution.  Results  are 
presented  for  constrained  and  unconstrained  water  volume  fraction  distributions  over  increasing 
levels  of  attenuation.  New  adjoint-based  algorithms  were  constructed  which  leave  the  final  time 
free  and  satisfy  all  first  order  necessary  conditions  as  well  as  avoid  taking  a  variation  at  the 
shock  front. 
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CHAPTER  1 : 
Introduction 


1.1  Ignition  Overpressure 

Ignition  overpressure  (IOP)  is  a  phenomenon  present  at  the  start  of  an  ignition  sequence  in 
launch  vehicles  using  solid-grain  propellants.  When  the  grain  is  ignited  the  pressure  inside  the 
combustion  chamber  quickly  rises  several  orders  of  magnitude.  This  drives  hot  combustion 
products  toward  the  nozzle  and  out  to  the  open  atmosphere  at  supersonic  speeds.  An  IOP  wave 
is  an  unsteady  shock  wave  which  originates  from  the  exit  plane  of  the  nozzle  and  propagates 
spherically  outward  near  sonic  conditions.  From  previous  launch  data  [l]-[3]  it  is  known  that 
overpressures  that  the  body  of  the  rocket  experiences  are  of  the  order  2:1,  or  about  15  psi.  The 
region  below  the  nozzle,  such  as  the  launch  pad  trench,  will  experience  further  compression  due 
to  displacement  of  gas  along  the  blast  wave’s  direction  of  propagation  and  overpressures  can  be 
as  high  as  10:1,  or  about  130  psi.  The  portions  of  the  IOP  wave  that  become  incident  on  the 
rocket  body  or  launch  platform  components  must  have  an  overpressure  below  a  known  thresh¬ 
old  to  avoid  costly  damage  and  enable  a  prolonged  lifetime.  This  is  the  purpose  of  the  water 
suppression  system  which  is  integrated  into  the  launch  platform.  The  current  technique  used 
by  NASA  and  other  launch  providers  is  to  spray  water  into  the  region  around  the  nozzle  before 
ignition.  This  forces  the  IOP  wave  to  propagate  through  water  before  becoming  incident  on 
the  rocket  body  or  platform  components.  Through  several  dissipative  mechanisms  this  causes  a 
sufficient  decrease  to  the  pressure  jump  across  the  shock,  preventing  damage. 

1.2  Experimental  Results  on  Droplet-Shock  Interactions 

A  survey  of  empirical  data  in  the  literature  was  carried  out  [4]-[21]  to  understand  how  droplets 
can  be  used  as  a  control  action  in  two-phase  flow.  With  blast  mitigation  as  the  process  of  interest, 
an  NRL  report  [7]  isolated  several  desirable  and  undesirable  effects.  Drag  on  the  droplets  by 
the  surrounding  gas  dissipates  momentum.  As  droplets  breakup,  the  gas  does  work  against 
the  surface  tension  which  dissipates  energy  from  the  gas.  Droplets  absorb  sensible,  latent  and 
radiative  heat  which  further  transfers  energy  from  the  gas.  Since  water  vapor  has  a  higher  heat 
capacity  than  air,  after  water  changes  phase  from  liquid  to  gas  the  resulting  gas  mixture  can 
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further  absorb  heat  from  the  gas.  An  undesirable  effect  occurs  when  liquid  turns  to  water  vapor 
so  rapidly  that  the  gas  mixture  density  increase  dominates  the  effect  on  pressure  according 
to  the  constitutive  relation.  The  relative  importance  of  each  interaction  mechanism  isn’t  fully 
understood  and  will  change  depending  on  the  flow  regimes. 

The  flow  regime  of  interest  will  have  shocks  present  with  pressure  ratios  (2:1-10:1)  and  tem¬ 
peratures  up  to  several  thousand  Kelvin.  Practical  droplet  sizes  have  a  diameter  Dp  less  than 
500p.m.  A  shock  incident  on  a  cloud  of  water  droplets  will  pass  through  them  so  quickly  that 
the  droplets  will  not  appear  to  react  until  after  the  shock  front  is  spatially  isolated  downstream. 
Droplets  will  start  to  get  dragged  along  in  the  direction  of  the  gas  flow,  breakup  and  vaporize. 
If  droplets  are  large  Dp  >  100pm,  then  in  the  presence  of  strong  shocks  they  catastrophically 
breakup  into  many  small  droplets  Dp  <  25pm.  Droplet  breakup  ceases  below  a  Weber  number 
of  12,  when  inertial  forces  dominate  over  surface  tension  [5]-[8].  Shock  attenuation  will  be 
greatest  when  the  droplets  can  sink  the  most  energy  out  of  the  gas  in  a  given  interval  of  space. 
This  will  cause  the  greatest  decrease  in  pressure  to  the  driving  gas  behind  the  shock.  Pressure 
information  will  travel  at  the  local  sound  speed  toward  the  shock  front  and  cause  a  decrease  in 
pressure  equal  to  the  diminished  driver  gas  pressure. 

There  are  general  trends  in  the  data  which  help  categorize  the  relative  importance  of  all  of  the 
dissipative  mechanisms  present.  First  and  foremost,  droplet  vaporization  can  extract  orders 
of  magnitude  more  energy  from  the  gas  as  droplet  breakup  can.  Secondly,  the  latent  heat  of 
vaporization  is  the  most  significant  dissipative  mechanism  [8].  Lastly,  it  has  been  shown  [12] 
that  a  non-dimensional  droplet-flow  parameter  can  predict  overpressure  attenuation  level,  shown 
in  Figure  1.1.  The  horizontal  axis  is  a  non-dimensional  flow  parameter.  The  vertical  axis 
is  the  decrease  in  overpressure  for  a  shock  interacting  with  droplets  divided  by  the  decrease 
in  overpressure  when  no  droplets  are  present,  as  measured  at  a  fixed  location  downstream. 
Shock  tube  experiments  were  conducted  which  varied  the  droplet  size  and  shock  strength.  The 
data  suggests  that  maximizing  the  exposed  surface  area  of  the  droplets  will  yield  the  greatest 
enhancement  of  IOP  attenuation. 

1.3  Recent  Research  on  IOP  Attenuation 

For  many  years,  the  water  injection  strategy  for  handling  the  shuttle’s  IOP  transient  blast  has 
been  based  on  order-of-magnitude  estimates  relating  the  amount  of  energy  dissipation  required 
to  the  total  amount  of  water  used  [22].  More  than  enough  water  was  used  and  the  results  were 
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Figure  1.1:  Experimental  data  on  IOP  attenuation,  Jourdan  et  al.  [12] 

sufficient.  For  the  heavy  lift  vehicles  of  the  future,  IOP  blast  waves  will  be  more  substantial 
and  require  a  better  understanding  of  how  the  water  affects  the  IOP  strength.  More  recently, 
CFD  research  is  underway  at  NASA  [23]-[26]  and  in  the  private  sector  to  predict,  with  greater 
precision,  the  launch  environment  during  ignition  for  various  solid  grain  rockets  [27,  28]. 

The  first  work  [29]  on  optimizing  one  of  these  precise  CFD  simulations  was  a  parametric  study 
that  looked  at  water  arrangement  in  the  nozzle  region  and  how  it  affected  the  maximum  IOP 
strength.  At  NASA  Huntsville,  Cannabal  showed  in  his  dissertation  [29]  and  a  later  publication 
[30]  that  attenuation  is  very  insensitive  to  droplet  velocity  and  demonstrated  the  existence  of  an 
optimal  water  injection  arrangement.  Water  cooling  the  plume  near  the  nozzle  has  the  greatest 
desirable  effect  of  attenuating  the  transmitted  IOP  strength.  However,  an  excessive  amount  of 
water  near  the  nozzle  causes  obstruction  to  the  blast  wave  and  intensifies  pressure.  The  results 
suggest  an  optimal  arrangement  of  water  exists  but  there  are  still  an  infinite  number  of  possible 
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water  distributions  even  in  one  dimension.  Trial  and  error,  or  cost  gradient  methods  based  on 
a  few  discrete  inputs  are  the  only  option  and  can  yield  only  coarse  notions  about  continuous 
optimal  water  distributions. 

1.4  New  Contributions  of  Current  Work 

The  objective  of  this  work  is  to  develop  a  computational  tool  that  can  directly  calculate  a  droplet 
optimal  control  for  attenuating  a  range  of  blast  waves  to  a  desired  minimal  overpressure.  As 
shown  in  Figure  1.2,  the  water  droplets’  control  action  will  be  the  result  of  the  initial  distribu¬ 
tion  of  the  water  volume  fraction  variable  oti(x,  0)  and  the  initial  droplet  diameter.  The  optimal 
control  a*(x,  0)  will  be  the  distribution  of  water  droplets  which  yields  the  greatest  decrease  to 
the  jump  in  overpressure  of  the  transmitted  shock  while  not  using  any  more  water  than  neces¬ 
sary.  Since  the  two-phase  system  and  resulting  control  problem  are  quite  complicated,  first  a 
single -phase  control  problem  is  formulated  where  the  control  action  is  a  distributed  energy  sink 
behind  the  shock,  a  simplification  for  the  way  droplet  vaporization,  the  dominant  dissipative 
mechanism,  affects  the  gas  flow.  A  more  sophisticated  two-phase  model  [31]  is  presented  and 
an  analogous  control  problem  is  formulated  in  which  the  control  takes  the  form  of  the  free  initial 
distribution  of  water  and  evolves  dynamically  according  to  the  two-phase  model.  These  control 
problems  are  new  formulations  and  the  algorithmic  solutions  developed  in  order  to  satisfy  all 
necessary  conditions  have  aspects  new  to  computational  results  for  unsteady  shock  wave  atten¬ 
uation.  Given  the  mathematical  framework  presented  in  Chapter  2  and  the  numerical  methods 
given  in  Chapter  3  of  this  dissertation,  engineering  results  such  as  those  presented  in  Chapter  4 
can  be  obtained  for  a  given  incident  blast  wave  boundary  condition. 
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Figure  1.2:  Physical  diagram  of  simulated  interaction. 
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CHAPTER  2: 

Mathematical  Theory 

2.1  Euler  System  in  One  Dimension 

The  mathematical  theory  of  hyperbolic  conservation  laws  has  been  well  established  in  the  20th 
century  [32]-[36].  In  order  for  a  PDE  to  be  an  accurate  physical  system  the  right  conserved 
quantities  must  be  identified  and  an  entropy  condition  must  exist.  The  conservation  equations 
dictate  that  the  rate  of  change  of  a  conserved  quantity  within  a  fixed  volume  equal  the  flux  of 
that  quantity  across  the  volume’s  surface.  The  entropy  condition  selects  which  conditions  are 
physical  when  discontinuous  solutions  are  present. 

2.1.1  Conservation  Equation 

The  time-dependent  Euler  Equations  are  a  system  of  non-linear  hyperbolic  conservation  laws 
that  describe  the  dynamics  of  compressible  fluids  where  the  effects  of  viscosity,  heat  conduction 
and  body  forces  are  negligible.  In  one  spatial  dimension  there  are  three  conserved  quantities: 
mass  density  p,  linear  momentum  pu  and  total  energy  pE  per  unit  volume. 

ltu+lF^  =  °  ai) 

U  =  (p(x,  t),pu(x ,  t),pE(x ,  t))T  (2.2) 


The  state  vector  of  conserved  quantities  U  is  a  function  of  a  time  variable  which  exists  in  the 
interval  [0,  T]  and  a  single  space  variable  x  which  exists  in  a  one  dimensional  domain  Q.  F(U) 
is  the  non-linear  flux  of  conserved  quantities. 


F(U) 


( 


pu 

pu 2  +  P 


\ 


V  u{pE  +  P)  ) 


(2.3) 
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For  an  ideal  gas,  the  constitutive  relation  that  relates  pressure  to  internal  energy  is  the  ideal  gas 
law. 


P  =  pe  (7  -  1) 


(2.4) 


c 

The  ratio  of  specific  heat  is  defined  7  =  ^  where  Cp  is  the  specific  heat  of  the  gas  at  constant 
pressure,  divided  by  that  at  constant  volume  Cv,  It  will  be  convenient  to  have  the  relation  be¬ 
tween  internal  energy  and  temperature  for  a  calorically  perfect  gas,  e  =  CVT.  A  final  definition 
relating  internal  and  total  energy  closes  the  system  from  which  all  thermodynamic  quantities 
can  be  calculated. 


pE 


(2.5) 


2.1.2  Quasilinear  Form 

Equation  2.1  can  be  written  in  quasi-linear  form  by  creation  of  the  Jacobian  matrix  A—Wj. 


dUi 

dt 


+  Aij(U) 


dUi 

dx 


0 


(2.6) 


The  elements  of  the  Jacobian  matrix  are  given  in  Equation  2.7. 


A(U) 


( 


0 


(pu)2  7  ~  3 

p2  2 

~7  (pu)  ( pE ) 

V  P2 


1  0  N 

-^(7-3)  (7-1) 

7  (pE)  3^  ^  (H2  7  {pu) 

— —-2(7-1)  -o-  — — 


(2.7) 
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2.1.3  Primitive  Form 


The  Euler  system  can  be  written  in  the  primitive  variable  basis  U  =  ( p,u,P )T  by  left  mul¬ 
tiplying  Equation  2.6  by  M-1  =  and  inserting  MM -1  between  the  Jacobian  matrix  A 
and  the  term  |^.  u  is  the  linear  gas  velocity  and  P  is  the  gas  pressure.  Defining  the  matrix 
A  =  AM1  AM  gives  the  primitive  form  of  the  Euler  system  in  Equation  2.8.  Details  of  the 
transformation  between  bases  for  a  single  dimensional  fluid  flow  are  given  in  Merkle  [39]. 


dUi 

dt 


+  A 


v 


0 


The  elements  of  A  are  given  in  Equation  2.9. 

0  ^ 
1/P 

u  ) 


A(u)  = 


u  p 
0  u 
^  0  yP 


(2.8) 


(2.9) 


2.1.4  Characteristic  Form 

It  can  be  shown  that  A  has  eigenvalues  (u,  u  +  c,  u  —  c)  where  c  =  y/'yRT  is  the  speed  of 
sound  in  the  gas.  Any  system  is  hyperbolic  if  the  eigenvalues  of  the  Jacobian  matrix  are  real 
and  distinct.  Clearly  this  is  true  in  the  case  of  the  Euler  System  in  ID.  The  specific  gas  constant 
of  air  is  denoted  as  R  and  the  gas  temperature  as  T.  When  the  corresponding  left  eigenvectors 
multiply  dU  the  characteristics  equations  are  obtained. 


dP  —  pcdu  =  0  along 
dP  —  c2dp  =  0  along 
dP  +  pcdu  =  0  along 


dx/dt  =  u  —  c 
dx/dt  =  u 
dx/dt  =  u  +  c 


(2.10) 


Solving  the  Riemann  Problem  for  the  unsteady  Euler  Equations  amounts  to  constructing  the 
rest  of  the  flow  field  between  these  characteristic  waves,  either  approximately  or  iteratively  up 
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to  arbitrary  precision  [40] . 


The  method  of  characteristics  takes  advantage  of  the  property  that  hyperbolic  systems  have 
distinct  finite  wave  speeds.  If  an  observer  is  traveling  at  a  characteristic  speed  the  observer  will 
measure  an  invariant  quantity  of  the  flow.  The  ID  Euler  system  in  characteristic  form  is  given 
in  Equation  2.11.  It  is  obtained  in  a  way  analogous  to  the  transformation  between  Equations 
2.6  and  2.8.  The  matrix  A  is  diagonal  and  it’s  elements  (eigenvalues)  Xt  are  the  real  and  distinct 
wave  speeds  (u,u  +  c,u  —  c). 


dU  AdU 
dt  dx 


0 


(2.11) 


Traveling  at  the  characteristic  speeds,  an  observer  would  measure  the  Riemann  invariants  U  = 
(S/R,Tc  +  u,Tc  —  u)T  as  constant.  S  is  the  entropy,  R  is  the  specific  gas  constant,  c  is  the 
speed  of  sound  and  T  =  2/ (7  —  1).  Traveling  at  the  speed  of  a  shock  wave  means  that  the  shock 
will  never  traverse  the  observer  and  hence  the  observer  will  never  see  an  increase  in  entropy. 
Indeed,  the  Riemann  invariants  are  useful  for  analyzing  isentropic  portions  of  compressible 
flows  such  as  rarefaction  waves. 


2.1.5  Well-Posedness 

For  hyperbolic  conservation  laws,  a  well-posed  problem  requires  real  and  distinct  eigenvalues 
of  the  system’s  Jacobian  matrix  [37].  This  is  equivalent  to  the  existence  of  a  symmetrizer  matrix 
S  that  has  real,  positive  definite  elements  and  has  the  property  that  SA  is  symmetric.  It  can  be 
shown  that  the  ID  Euler  System  has  a  symmetrizer  matrix  as  long  as  density,  or  temperature 
remain  positive  [43].  This  assumption  is  built  in  to  the  notions  of  either  quantities  therefore  the 
system  is  always  symmetrizeable  and  hence  the  Riemann  Problem  for  the  Euler  Equations  is 
well-posed. 

2.1.6  Entropy  and  Shocks 

From  the  definition  of  entropy  and  the  primitive  form  of  the  mass  and  momentum  conservation 
equations  it  can  be  shown  that  entropy  is  constant  along  dx/dt  =  u.  That  is,  regions  of  smooth 
flow  are  isentropic.  The  presence  of  a  discontinuity  in  velocity  necessitates  that  a  shock  wave 
has  formed  where  the  jump  in  flow  quantities  are  determined  uniquely  by  the  Rankine-Hugoniot 
Jump  Conditions.  For  a  shock  wave  moving  left  to  right,  denote  the  state  just  behind  the  shock 
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Ul  and  the  state  immediately  upstream  of  the  shock  Ur.  The  jump  quantities,  denoted  with 
square  brackets  [  ],  are  given  in  Equation  2.12. 


OF  (U)  [F  ([/)]  =  F  (UR)  -  F  (UL) 
d(U)^[U]  =  UR-UL 


Rearranging  Equation  2.12  gives  the  shock  speed  in  terms  of  the  jump  across  the  shock.  Equa¬ 
tion  2.13  is  the  Rankine-Hugoniot  jump  conditions  obtained  without  any  specification  to  a  fluid. 
This  condition  is  a  property  of  any  shock  wave  in  a  hyperbolic  system. 


dx 

~ 77  shock 


[F  m 
[u] 


(2.13) 


Applying  Equation  2.13  to  the  Euler  Equations  moving  at  the  shock  speed  gives  the  jump  in 
primitive  variables  in  a  coupled  system  of  algebraic  equations  which  can  be  solved  iteratively 
as  will  be  described  in  Chapter  3.  In  addition,  an  alternative  numerical  method  of  Godunov  will 
be  described  which  is  conservative  and,  with  a  stability  criteria,  obeys  Equation  2.13  without  an 
iterative  solver. 

Since  the  eigenvalues  are  increasing  functions  of  the  state  variables,  the  the  flux  is  convex  with 
respect  to  U.  Lax  [32], [33]  showed  that  this  requires  characteristics  on  either  side  of  a  shock  to 
run  into  that  shock.  For  a  system  with  a  right  moving  shock  the  entropy  condition  is  given  in 
Equation  2.14  for  k  number  of  waves  where  1  <  k  <  n,  and  n  is  the  number  of  characteristic 
waves. 


dx 

^k—  liJJ l)  ^  ~q^\ shock  ^  Lj) 

dx 

\k(UR)  <  —\ shock  <  ^k+l(UR) 


(2.14) 


Intuitively,  this  agrees  with  the  physical  notion  that  entropy  only  increases.  For  the  ID  Euler 
system  n  =  3.  Since  a  shock  requires  supersonic  conditions  the  shock  speed  must  be  faster 
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than  the  fastest  rate  at  which  any  information  can  travel  in  the  ambient  upstream  Right  state, 
so  it  must  be  true  that  k  =  3.  This  means  that  the  entropy  condition  also  requires  exactly  one 
characteristic  in  the  Left  state  behind  the  shock  to  be  faster  than  the  shock  itself.  The  shock 
is  driven  by  the  upstream  pressurized  gas  which  must  communicate  its  presence  via  the  u  +  c 
characteristic  in  the  Left  state.  This  is  not  hard  to  believe  since  the  jump  conditions  show  that 
static  temperature  increases  across  the  shock  which  raises  the  speed  of  sound  in  the  Left  state. 
Figure  2.1  shows  the  instantaneous  wave  pattern  of  a  shock  wave  formed  from  a  discontinuous 
and  stationary  initial  condition. 


Figure  2.1:  Instantaneous  wave  patterns  of  a  Left  and  Right  state  and  shock  wave  which  obey 
the  entropy  condition  of  the  Euler  system  in  ID 

From  the  first  law  of  thermodynamics  the  entropy  between  states  ’L’  and  ’R’  is  given  in  Equation 
2.15. 


Sl  —  Sr  =  —R  ■  In 


Pol 

Pqr 


(2.15) 


Pol  and  P0R  are  the  total  pressures  to  the  left  and  right  of  the  shock  respectively.  From  manip- 
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ulation  of  Rankine  Hugoniot  Jump  conditions  it  can  be  shown  that  for  supersonic  flows  there  is 
necessarily  a  decrease  in  total  pressure,  verifying  by  Equation  2.15  that  indeed.  Lax’s  mathe¬ 
matic  entropy  condition  on  the  Euler  system  in  one  spatial  dimension  yields  the  thermodynamic 
principle  that  entropy  can  only  increase  across  a  shock. 

2.1.7  The  Riemann  Problem 

The  Riemann  Problem  is  an  initial  value  problem  where  the  initial  data  has  a  discontinuity. 


U{x,Q) 


UL  x  <  0 
Ur  X  >  0 


(2.16) 


For  gas  dynamics,  the  Riemann  Problem  is  the  general  form  of  the  shock  tube  problem  [40]. 
Two  stationary  gases,  one  at  significantly  higher  pressure  and  possibly  temperature,  are  sepa¬ 
rated  by  a  thin  wall.  As  the  thin  wall  ruptures,  a  shock  wave  propagates  into  the  gas  held  at 
lower  pressure  while  a  rarefaction  wave  travels  through  the  gas  held  at  higher  pressure.  Initially 
a  density  discontinuity  exists  between  the  two  gases  and  this  contact  surface  propagates  into 
the  lower  pressure  gas  at  the  speed  determined  by  the  jump  in  velocity  at  the  shock.  The  shock 
travels  at  the  speed  of  the  gas  behind  it  plus  the  speed  of  sound  of  the  low  pressure  gas.  The 
rarefaction  wave  travels  in  the  opposite  direction  of  the  gas  at  a  speed  equal  to  the  speed  of 
the  gas  minus  the  speed  of  sound  in  the  high  pressure  gas.  The  wave  pattern  of  the  shock  tube 
problem  at  the  instant  the  diaphragm  bursts  is  shown  in  Figure  2.2. 

It  is  worth  pointing  out  that  the  fastest  wave  speed  may  be  u  —  c  since  the  speed  of  sound  may  be 
much  greater  in  the  high  pressure  gas.  In  Chapter  3,  the  maximum  wave  speed  will  determine 
the  maximum  allowable  timestep  discretization  and,  therefore,  computing  expense. 

By  algebraic  manipulation  of  the  Rankine-Hugoniot  Jump  Condition  and  the  isentropic  equation 
of  state  for  an  ideal  gas  it  can  be  proven  [40]  that  solving  the  Riemann  Problem  for  the  Euler 
system  amounts  to  finding  the  rootp*  of  Equation  2.17. 


f(p*,  Ul,  Ur)  =  fi (p* ,  Ul)  +  fn(p*,  Ur)  +  [U]  =  0  (2.17) 

The  rest  of  the  primitive  variables  u*,p*  can  be  uniquely  calculated  from  p*  and  the  correct 
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Star  Region 


Figure  2.2:  Solution  regions  of  the  shock  tube  problem  and  wave  pattern  at  the  initial  instant 
the  diaphragm  bursts 


definition  of  Jl  and  Jr.  The  present  analysis  will  be  restricted  to  a  left  moving  rarefaction  wave 
and  a  right  moving  shock  wave.  The  isentropic  relation  between  pressure  and  density  holds 
across  rarefaction  waves.  Solution  to  the  flow  in  the  U*l  region  is  derived  from  the  constant 
Ul  state.  By  combining  the  isentropic  relation  and  the  Riemann  invariants  fjJ  for  a  left  moving 
rarefaction  can  be  written  in  terms  of  the  constant  left  state  Ul  and  the  free  variable  p*  in 
Equation  2.16.  cl  is  the  speed  of  sound  in  the  constant  U l  state. 


=  :)  ai8) 

f[{  in  the  U*r  region  is  determined  by  the  constant  Ur  state  and  the  Rankine-Hugoniot  jump 
conditions.  Applying  Equation  2.13  to  the  ID  Euler  Equations  in  a  frame  of  reference  moving 
along  with  the  shock  gives  the  jump  conditions  and,  with  some  rearranging,  yields  fR  for  a  right 
moving  shock  wave. 
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1/2 


(2.19) 


h(p,UR)  =  (f-pR )  (j+-) 

The  constants  Ar  and  Br  are  given  in  Equation  2.20.  The  flow  can  be  uniquely  determined  at 
any  time  in  the  star  region  given  a  physical  initial  condition. 


Ar  = 


2 

(7  +  1)  PR 


Br  = 


(7-1) 
(7  +  1) 


PR 


(2.20) 


The  exact  solution  requires  an  iterative  procedure  at  every  sample  point  in  space  and  at  any 
desired  moment  in  time.  This  is  not  practical  for  large  domains  or  for  the  added  complexity  of 
another  phase.  The  difficulty  in  accuracy  of  any  approximate  solution  to  the  Riemann  Problem 
comes  from  any  discontinuities  at  shock  fronts.  Here  finite  difference  methods  fail  and  a  finite 
volume  method  based  on  the  integral  form  of  the  conservation  equation,  rewritten  in  Equation 
2.21,  is  better  suited. 


U (x,  t2)dx 


rx  2  rt2  rt2 

/  U(x,ti)dx  +  /  F(U(xi,  t))dt  —  /  F(U(x2,t))dt 

lx  1  Jt\  Jt\ 


(2.21) 


No  discretization  or  approximation  has  been  made  in  Equation  2.21  and  it  is  true  for  any  general 
domain  [x1,x2\  x  [ti,  t2]  - 

Numeric  schemes  based  on  primitive  variables  fail  at  shock  waves.  If  the  jump  in  quantities 
across  the  shock  are  wrong,  the  shock  speed  will  be  wrong  and  therefore,  over  time,  the  lo¬ 
cation  will  be  wrong  as  well.  Only  schemes  based  on  conservative  variables  have  acceptable 
accuracy  near  shocks.  A  numeric  method  is  conservative  only  if  the  approximate  numerical  flux 
is  constructed  in  such  a  way  that  each  quantity  is  conserved  at  every  time  step.  In  particular, 
Lax  [36]  showed  that  conservative  methods  that  are  also  convergent  will  converge  to  the  cor¬ 
rect  weak  solution  at  shocks.  With  that  in  mind,  Godunov’s  Method  with  a  suitable  choice  of 
numerical  flux  will  be  described  in  Chapter  3  and  used  for  all  results  shown  in  Chapter  4. 
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2.2  Two-Phase  Flow  Model 


2.2.1  Assumptions 

The  vast  majority  of  the  two-phase  flow  model  is  based  on  the  work  of  Saurel  and  Abgrall 
[31], [41].  It  consists  of  six  balance  laws  and  a  seventh  non-conservative  PDE  for  the  volume 
fraction  variable.  The  model  was  chosen  because  it  maintains  a  hyperbolic  structure  and  is  ap¬ 
plicable  to  physical  phenomena  of  interest,  water  droplet-shock  interactions.  This  is  achieved 
by  considering  both  phases  or  fluids  as  compressible,  by  assuming  many  droplets  per  cell  and 
therefore  homogenized  interface  conditions  and  by  solving  the  same  set  of  equations  every¬ 
where  in  space. 

2.2.2  Balance  Equations 

The  conservative  vector  of  the  flow  is  again  denoted  as  U.  The  subscripts  in  Equation  2.22 
denote  gas  or  liquid. 


U  {Oig,  OigPg,  Ol  g P  gU  g  ,  Ol  g  g  E g  ,  Ol  /  /  ,  CpP/M/,  Cp  P/ E)  ) 


(2.22) 


dOLg 

~dt 


+  Vt 


OOLg 

dx 


0 


(2.23) 


d_ 

Ft 


/ 


V 


OLgPg  ^ 

agPgug  ^ 

(  o  N 

agPgug 

agPgU2g  +  agPg 

Pi 

agPgPg 

d 

+ 

ug  ( agPgPg  +  agPg ) 

PrVi 

dag 

Oil  pi 

dx 

oiipm 

0 

dx 

onpm 

OilPiu 1  +  aiPi 

~Pi 

onpiEi  j 

y  ui  ( aipiEi  +  aiPi )  ) 

y  -PF  ) 

+  S{U) 


(2.24) 


The  gas  phase  obeys  the  ideal  gas  equation  of  state  while  the  stiffened  gas  equation  of  state  is 
used  for  the  liquid. 


Pgeg 


Pg 

7-1 


(2.25) 
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Pi  +  ng  i 

li  -  1 


(2.26) 


Piei  = 


7  =  1.4  is  the  gas  constant  of  air,  7/  =  4.4  is  the  analogous  constant  for  water  [31]  and 
77  =  6  •  108  Pa  is  the  stiffening  constant  that  makes  large  changes  in  liquid  pressure  produce 
almost  no  changes  in  density.  Equations  2.27  and  2.28  are  closure  relations  for  the  internal 
energy  pe  to  the  total  energy  pE,  per  unit  volume,  for  each  phase. 


PgEg  PgCg  +  ^  PgUg 


(2.27) 


Pi  Ei  =  piei  +  -piuf 


(2.28) 


2.2.3  Interface  Quantities 

The  volume  fraction  will  propagate  at  a  mean  inter-facial  velocity  which  is  a  center-of-mass 
estimate  given  in  Equation  2.29. 


y -  _  agPgug  P  aiPiui  ^  29) 

agPg  T  OliPi 

Ei  and  P,  are  volume  averages  of  total  energy  and  pressure  respectively  at  the  interface  shown 
in  Equation  2.30. 


Ei  otgEg  T  aiE[ 

Pi  agPg  ^ iPl 


(2.30) 


The  sum  of  the  volume  fractions  of  each  phase  will  always  equal  1  for  all  x  £  Q. 


OLg  T  Oil  -  1 


(2.31) 


When  either  volume  fraction  tends  toward  zero,  the  ID  Euler  system  is  recovered. 
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2.2.4  Primitive  Form 


As  with  the  single-phase  Euler  system,  the  two-phase  system  can  be  written  in  non-conservative 
form.  Again  let  U  be  the  primitive  vector. 


U  (Q!,g:  Pg,  Ug,  Pg,  Pl,  Ul, 


(2.32) 


In  this  basis  the  Jacobian  matrix  for  the  system  of  Equations  2.23  and  2.24  is  given  in  Equation 
2.33. 


A(U)  = 


V 


Vi 

0 

0 

0 

0 

0 

0 

Pg_ 

«9 

{V- 

U9 ) 

Ug 

Pg 

0 

0 

0 

0 

Pg-P. 

agPg 

0 

Ug 

1/Pg 

0 

0 

0 

Pscl 

ag 

HVi- 

~U9 ) 
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Pgpg 

Ug 
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0 

0 

PL 

Oil 

(Vi- 

Ul ) 

0 

0 

0 

Ul 

pi 

0 

9 

Pi-Pi 

OigPg 

0 

0 

0 

0 

Ul 

1  /pi 

pK 

‘(Vi- 

-  Ul) 

0 

0 

0 

0 

pi  A 

Ul 

J 


(2.33) 


The  source  term  H (U)  which  multiplies  in  Equation  2.24  has  been  coupled  into  the  Jaco¬ 
bian  matrix  and  the  two-phase  quasi-linear  system  now  has  the  form  shown  in  Equation  2.34. 


f  +  MU)aP  =  Sl{u) 


(2.34) 


With  this  fortunate  manipulation,  the  seven  eigenvalues  can  be  found  by  diagonalizing  A(U). 
The  characteristic  speeds  are  (V),  ug,ug  —  cg,ug  +  cg,ui,ui  —  ci,ui  +  q),  where  q  is  the  local 
speed  of  sound  in  the  liquid  phase.  The  velocity  of  the  phase  interface  is  the  seventh  charac¬ 
teristic  wave  speed  of  the  system.  All  are  distinct  except  locally  where  they  are  degenerately 
zero. 


2.2.5  Phase  Interactions 

The  source  vector  S(U)  is  broken  up  into  three  separate  interactions  as  defined  in  Equation 
2.35. 
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S(U)  =  MH(U)  +  VR(U)  +  PR(U) 


{  0  1 

f  0  l 

1  ,i(p,-Pi)  ) 

m 

0 

0 

mVi 

Fd 

0 

?Tl  (EhV  +  Pi)  +  Qi 

+ 

PdVi 

+ 

-pPi  (Pg  -  Pi) 

—tin 

0 

0 

—mVi 

~Fd 

0 

y  di'  (Pfiv  T  Pi)  Qi  J 

v  -FdVi  ) 

v  R pi  ( Pg  -  Pi)  ) 

(2.35) 


The  source  vector  for  mass  and  heat  exchange  between  the  phases  is  denoted  as  MH(U ).  The 
source  terms  for  velocity  and  pressure  equilibration  are  VR(U )  and  PR(U )  respectively.  As 
Saurel  and  Abgrall  point  out  [31],  an  interface  separating  two  phases  must  reach  the  same 
pressure  through  microscopic  interactions.  Indeed,  without  enforcing  this  condition,  notions  of 
thermodynamic  properties  such  as  temperature  cannot  be  determined  and  numerical  oscillations 
due  to  pressure  differences  will  grow  without  significant  artificial  dissipation. 

The  equilibrium  condition  Pg  =  I)  is  chosen  thereby  neglecting  the  effect  of  surface  tension. 
As  stated  in  Chapter  1,  the  effect  of  droplet  breakup  is  omitted  in  favor  of  vaporization  since  the 
latter  can  be  shown  to  be  much  more  significant  of  an  energy  sink  to  the  gas.  The  microscopic 
pressure  equilibration  causes  a  volume  and  internal  energy  variation  of  each  phase.  Details  of 
the  solution  procedure  for  the  two-phase  model  are  saved  for  Chapter  3.  Isolating  the  ODE 
=  PR(U)  in  Equation  2.35  yields  conditions  better  suited  to  computation. 


OOig 

dt 


P,  (. Pg  ~  Pi) 


(2.36) 


dOLgpgEg 

dt 


^ R  ( Pg  ~  Pi) 


(2.37) 


da  g  pt  Ei 
dt 


~pPi  (Pg  -  Pi) 


Therefore,  Equations  2.36-2.38  simplify  to  Equations  2.39  and  2.40 


(2.38) 
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(2.39) 


DOLgpgEg  _  'dOlg 

dt  1  dt 


dOigplEl  _  d  Qtg 

dt  ~  dt 


(2.40) 


In  a  closed  system,  this  is  just  the  first  law  of  thermodynamics  for  an  isentropic  transformation. 
Time  integration  of  both  sides  yields  conditions  given  in  Equations  2.41  and  2.42  whose  solution 
will  be  discussed  in  Chapter  3. 


(agpgEg)PR  -  ( agpgEg )°  =  f  9  Pidag  (2.41) 

Ja° 

(alPlEi)PR  -  {cnptEif  =  -  9  Pidoig  (2.42) 

The  superscript  PB  denotes  the  solution  to  the  ODE  for  pressure  relaxation  while  0  denotes  the 
starting  value  which  comes  from  the  solution  to  the  balance  PDE  at  the  current  time  step.  The 
drag  force  Fd  exerted  by  the  gas  onto  a  spherical  water  droplet,  which  is  responsible  for  the 
VFL{U )  source  term,  is  given  by  the  empirical  drag  law  in  Equation  2.43 

Fd  =  CdpgD2pai  ( ug  -  uif  (2.43) 

where  Dp  is  the  diameter  of  the  water  droplet.  This  is  a  dynamic  variable  distributed  in  space 
with  the  assumption  of  locally  mono-dispersed  droplets.  Cd  is  a  constant  drag  coefficient. 

The  exchange  of  mass  and  heat  between  the  phases  will  be  due  to  vaporization  of  liquid  water 
droplets  by  the  surrounding  gas.  The  rate  of  gaseous  mass  production  in  the  form  of  water  vapor 
from  the  liquid  phase,  m,  is  defined  in  Equation  2.44. 


rii  —  Ph2Ov  ( Pg )  di 


(2.44) 


The  rate  of  change  of  the  volume  fraction  of  water  within  a  fixed  volume  of  space  Ax3  is  related 
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to  the  diameter  rate  of  change  via  Equation  2.45. 


SP  dDp 
Ax3  dt 


(2.45) 


Sp  denotes  the  total  surface  area  of  droplets  per  a  fixed  volume  defined  in  Equation  2.46. 


S'p  =  47r(:y)  np 


(2.46) 


Np  is  the  number  of  droplets  per  the  same  fixed  volume  as  is  uniquely  determined  by  the  volume 
fraction  of  water  cq  for  mono-dispersed  droplets. 


Np  cq 


Ax3 


(2.47) 


The  evolution  of  droplet  size  needed  in  Equation  2.47  is  determined  by  the  rate  of  vaporization. 
In  this  work,  the  Empirical-Beta  Vaporization  law  is  implemented  with  constant  values  coming 
from  [7]. 


d  d2 

-gf  =  -P(Tg)  (2.48) 

(3  (Tg)  =  7600  (l  +  7.4  •  10“7  (Tg  -  300A')+)  //- m2 / s  (2.49) 

As  will  be  demonstrated  in  Chapter  4,  the  range  of  pressure  magnitude  for  the  flows  of  interest 
to  the  IOP  attenuation  problem  are  from  1-10  atm.  According  to  steam  tables  [42],  the  water 
vapor  density  will  vary  by  an  order  of  magnitude  over  this  pressure  range.  Consequently,  water 
droplet  vaporization  will  be  much  more  effective  to  high  pressure  shocks.  Equation  2.50  is  the 
least-squares  quadratic  fit  to  the  steam  table  data  for  water  vapor  density  as  a  function  of  the 
surrounding  gas  pressure. 


pH2Ov{Pg/atm )  «  -.0024344  •  P°g  +  .5368  •  Pg  -  .077246  kg/m 3  (2.50) 
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The  latent  heat  of  vaporization  of  water  at  its  boiling  point  Ljlv  =  2.26  MJ/kg. 

The  last  term  in  the  MH(U)  vector,  0,  is  the  rate  of  heat  exchange  between  phases  at  the 
interface  given  in  Equation  2.51. 


Qi  —  hSp  (Tg  —  Tt)  (2.51) 

h  is  the  heat  transfer  coefficient.  For  water  droplets  of  diameter  Dp,  h  =  — — —  where  Nu  is  the 

Up 

Nusselts  Number  and  A  =  0.6  W/mK,  is  the  thermal  conductivity  of  water. 

2.2.6  Well-Posedness 

The  Jacobian  matrix  in  Equation  2.33  is  symmetrizeable  if,  in  addition  to  the  restriction  in  the 
ID  Euler  System,  the  volume  fraction  of  each  phase  remain  non- zero  [31].  Since  our  model 
dictates  that  the  same  equations  be  solved  in  each  cell,  this  restriction  is  not  a  problem. 

The  eigenvalues  of  A  above  are  (Vi,  ug,  ug  —  cg,  ug  +  cg,  ui,  ui  —  ci,  ui  +  q).  From  the  defini¬ 
tion  for  interface  velocity  and  sound  speeds  it  can  be  seen  that  all  these  eigenvalues  are  real 
and  distinct  if  a  few  conditions  hold.  If  there  are  regions  in  our  domain  where  a  single  phase 
is  completely  absent,  then  the  interface  velocity  will  be  degenerate  to  the  velocity  of  the  phase 
100%  present.  Therefore,  our  calculation  restricts  the  volume  fraction  of  either  phase  to  re¬ 
main  above  a  certain  small  threshold.  As  with  the  single-phase  Euler  system  in  one  dimension, 
well-posedness  can  be  shown  by  the  existence  of  a  real,  positive  definite  symmetrizer.  It  was 
shown  by  Texier  [43]  that  a  necessary  condition  on  the  positive-definiteness  of  the  symmetrizer 
depends  on  the  positiveness  of  a  certain  constitutive  variable  that  should  always  physically  be 
positive  such  as  density,  temperature  etc.  The  multi-phase  system  has  the  same  requirement.  In 
addition,  compressible  flow  has  been  shown  to  change  type  from  hyperbolic  to  elliptic  where 
the  eigenvalues  of  the  system  are  complex.  A  famous  example  is  the  sonic  bubble  in  aerody¬ 
namics  [44].  In  that  case,  the  flow  is  irrotational  and  allows  a  velocity  potential  formulation 
which,  in  polar  coordinates,  will  change  type  at  the  sonic  line.  Since  the  flow  is  steady,  a  steady 
sonic  bubble  develops  inside  in  which  the  flow  is  elliptic  and  hence  hyperbolic  flux-based  con¬ 
servative  numeric  methods  would  fail.  Blast  waves,  on  the  other  hand,  are  unsteady  and  not 
irrotational.  Any  sonic  regions  of  the  flow  are  dynamic  and  do  not  form  sonic  bubbles. 
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2.3  Optimal  Control  of  Systems  of  Partial  Differential  Equa¬ 
tions 

Optimal  control  theory  and  the  Pontryagin  Minimum  Principle  are  well-established  in  mathe¬ 
matical  and  applied  engineering  literature  [45] -[5  2]. 

2.3.1  Derivation  of  Necessary  Conditions  of  General  Optimal  System 

This  section  contains  a  general  derivation  for  all  the  necessary  first-order  conditions  of  opti¬ 
mality  when  the  dynamical  constraints  are  a  system  of  partial  differential  equations  with  free 
or  fixed  initial  and  final  data,  free  or  fixed  spatial  boundary  data  and  free  final  time.  In  Section 

2.4.2  the  general  conditions  are  applied  to  the  unsteady  Euler  Equations  in  ID  to  formulate  a 
shock  attenuation  problem  and  to  a  two-phase  flow  problem  in  Section  2.4.3  governed  by  the 
model  given  in  Section  2.2. 

To  start,  a  hyperbolic  system  of  balance  laws  under  a  distributed  control  action  zt  written  in  a 
general  conservative-flux  form  using  tensors  is  shown  in  Equation  2.52. 

a-w+kFhu)=sm+Zi  <2-52) 

where  i  =  1  where  m  is  the  number  of  state  variables  and  k  =  1, ....,  n  where  n  is  the 
number  of  spatial  dimensions.  If  the  derivatives  of  the  flux  vector  is  taken  with  respect  to  the 
state  vector  U,  a  Jacobian  matrix  Ak  for  each  of  the  klh  spatial  dimensions  is  defined. 

4(f/)  =  w,F'(u)  <2-53) 

The  state  dynamics  can  then  be  written  in  matrix  Jacobian  form.  It  is  important  to  mention  here 
that  writing  the  PDE  system  in  Jacobian  form  does  not  imply  that  this  is  the  form  used  to  solve 
the  flow.  For  flows  with  shock  waves,  conservative  methods  are  necessary  and  lend  themselves 
to  the  conservative- flux  form  of  the  equations.  However,  the  adjoint  system  of  equations  derived 
below  are  derived  from  the  Jacobian  form  and  are  always  linear  in  the  adjoint  variables,  making 
their  numerical  solution  easier. 
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Si{U)  +  Zi 


(2.54) 


dUi 

dt 


+  ^(U) 


9U± 

dxk 


A  non-linear  partial  differential  operator  N  can  be  defined  so  that: 


N(U,  z)  =  +  4 (U)^-  -  Si(U)  -Zi  =  0  (2.55) 

After  the  dynamical  system  is  known,  a  cost,  or  objective,  functional  must  be  intelligently 
defined.  In  the  most  general  framework  there  can  be  a  cost  associated  with  the  initial  data  I,  a 
cost  involving  the  final  data  K  and  a  running  cost  L  that  accumulates  over  the  time  interval  of 
the  control  problem  [0,  T] .  So  far,  initial  and  final  data  can  still  be  free  or  fixed,  no  restriction  has 
been  made.  The  cost  functionals  describe  the  objective  associated  with  how  these  parameters 
are  distributed.  Let  J  be  the  total  cost  functional. 


J  =  /([/(•,  0))  +  K(U(-,  T))  +  [T  L(U,  z)dt  (2.56) 

Jo 

To  incorporate  the  constraint  in  the  minimization,  N  is  multiplied  by  a  generic  continuous  vector 
V (x,  t )  and  added  to  J.  The  problem  then  becomes  minimizing  the  augmented  cost  functional 
J. 


J  =  /([/(•,  0))  +  K(U(-,  T))  +  [T  L(U,  z)  +  (V,  N(U,  z))  dt  (2.57) 

Jo 

It  is  standard  to  refer  to  Vi(x,t )  as  the  Lagrange  multiplier  or  adjoint  vector.  It  is  a  Tri¬ 
dimensional  vector  since  there  is  an  adjoint  state  for  every  state  variable. 

Let  (z*,  U*,V*,T*)  be  the  optimal  control,  state  vector,  adjoint  vector  and  final  time  respec¬ 
tively.  Let  a  perturbed  control  function  8z  be  added  to  the  optimal  control  solution  so  that: 


z  =  z*  +  e5z  (2.58) 

where  e>  0  is  a  small  constant.  The  perturbed  optimal  control  introduces  a  perturbation  to  the 
optimal  state  and  adjoint  state: 
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u  =U*  +  eSU 


(2.59) 


V  =  V*  +  eSV  (2.60) 

Another  way  to  see  5U  which  will  be  useful  shortly  is: 

SU  =  ^-U  (z* +  e8z)\e=0  (2.61) 

It  has  been  pointed  out  in  the  literature  [53] -[55]  that  if  the  state  variables  have  shocks,  a  pertur¬ 
bation  is  not  small  in  the  neighborhood  of  the  shock  and  does  not  have  the  vanishing  properties 
as  e  — »  0.  A  slight  increase  in  the  amplitude  behind  the  shock  perturbs  the  speed  shock  and 
therefore  also  the  location  of  the  shock  front.  This  causes  small  perturbations  to  induce  varia¬ 
tions  on  the  order  of  the  jump  across  the  shock.  The  presented  method  of  solution  avoids  this 
issue,  as  will  be  demonstrated  later.  Since  only  decreasing  the  amplitude  of  a  shock  wave  is 
desired  it  is  apparent  that  any  realistic  control  action  will  only  slow  the  shock  wave  down.  The 
target  state  Q(x)  and  final  time  penalty  will  be  constructed  in  such  a  way  that  all  variations 
of  the  solution  will  occur  upstream  of  the  shock  front  only.  Matching  the  simulated  pressure 
profile  under  control  action  with  the  target  final  state  near  the  shock  front  will  occur  by  allowing 
the  final  time  to  be  free.  Henceforth,  it  can  be  assumed  that  all  variations  are  taken  in  smooth 
regions  of  the  flow  and  that  the  solution  procedure  will  not  depend  on  a  shock  location  variable 
and  adjoint  state,  nor  is  a  more  sophisticated  variation  required. 

The  first  order  necessary  optimal  condition  that  must  be  satisfied  by  perturbed  control,  state, 
adjoint  state  and  final  time  (Sz,  SU,  SV,  St)  is: 

^-J(z*  +  e8z)  |e=0  =  0  (2.62) 
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=  £l  (U  (, z *  +  edz ))  |e=0  +  fj<  (i u  (z*  +  eSz))  |e=0 

So  L  (U  ( z *  +  edz) ,  z*  +  efc)  |e=o  +  (V*  +  edV,  N  (U  ( z *  +  edz) ,  z*  +  edz))  |e=o dt 

=  £l  (U  (z*  +  eSz))  |6=0  +  £K  (U  (z*  +  eSz))  |e=0 

+  £ L  (U  ( z *  +  edz) ,  z*  +  edz)  |e=o  +  £  ( V *  +  edV,  N  (U  (z*  +  edz) ,  z*  +  efc))  |e=odf 

+  (L+(y%7V(C/*,2*)))|t=T<5f  =  0 

=  (§,6U(x,oj)  +  (§£,  fl7(s,T))  +f0T  (§§,dU(x,t)) 

+  (§,  fc(a;,  t))  +(V*,N  ([/*,  1/*))  +  (V*  +  e<5V,  £ N  (U  ( z *  +  edz) ,  z*  +  edz))  \e=0dt 
+  (L+(V*,N(U*,z*)))\t=Tdt 


(The  second  term  in  the  second  line  is  zero  since  N  =  0). 

=  (§,5U(x,0j)  +  (%,6U(x,Tj) 

+  /oT  ( §U’^U(x,t ))  +  (^,Sz(x,t))  +  (V*,  £N  (U  (z*  +  edz) ,  z*  +  edz)  |e=0)  dt 
+  (L  +  (V*,N  (U*,  z*)))  \t=Tdt  =  0 


(2.63) 

This  is  a  good  place  to  see  structure  of  the  equations.  In  Equation  2.63  it  is  easier  to  see  the 
Frechet  derivative  analogously  as  a  regular  chain  rule  derivative.  It  is  worth  noting  that  in  this 
form  the  equation  is  valid  for  any  differential  operator  N,  linear  or  non-linear,  PDE  or  ODE. 
Note  that  the  components  of  the  objectives  I,  K,  L  are  clearly  separated  and  what  remains  is  to 
carry  the  perturbation  and  differentiation  with  respect  to  e  through  on  the  differential  operator 
and  then  integrate  by  parts  to  get  the  adjoint  partial  differential  system  with  boundary  condi- 
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tions.  That  term  in  the  above  equation  can  now  be  written  in  the  Jacobian  matrix  form  with  an 
assumption  that  the  eigenvalues  are  real  and  distinct. 


(V*,  fN  (U  (z*  +  eSz) ,  ^  +  e5z)  |e=0)  = 


{V*,  +A(V  +  e5z))  8U(z*+e8z)  +  efc))  _  (z*  +  eSz^  |£=o) 

=  (v,  a~w)  +  {v,  iA  (u  (z*  +  eSz))  du%yz))  1 6=0  ~  (v*,  §  su)  -  (V,6z) 


(2.64) 


The  last  line  shows  the  differential  terms  separated  and  makes  use  of  the  definition  of  SU.  The 
last  two  terms  in  Equation  2.64  are  simple  to  handle  but  the  first  two  require  more  work.  It  is 
now  necessary  to  use  tensor  notation. 


Term  1 


(l/-*  d  dUi(z*+e8z)\  \ 

[Vi  ’  de  di  )  le=0 


( t /*  8  dUj(z*+eSz)\  | 

\Vi  ’  at.  de  )  le=0 

(*?, 


dSUj 

dt 


=  [v:ml  -  (¥.«) 


(2.65) 
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Term  2 


(V*,  £Atj  (U  (z*  +  e8z )) 


dx  ) 


e=0 


{V",  wiiui  (**  +  ■*)  US'  +  4*  (C7  (**  +  £fc))  ilk  ("•  + £fe)  U 

{wtltV^SU,)  +  £  -  {^,SUs)  -  (Ay5,«/: 


(2.66) 


((wfS  -  +  £  (-W.Wi)  -  K??.^ 


Term  3 

ss- (c/  <*'• + £fc)>)  i-  =  wf°)  =  (ft^  (2'67) 

The  fourth  term  simply  yields  an  inner  product  between  the  adjoint  vector  and  5z.  It  will  make 
things  clearer  to  see  the  two  boundary  terms  in  isolation  from  the  differential  relations  [V*SU] q 
md£[(AV*,SU)}. 

After  integrating  by  parts  all  the  terms  in  the  last  inner  product.  Equation  2.63  can  be  written, 
for  hyperbolic  PDEs  with  source  terms,  as  [56]: 


(g»6U (x,0))  +  (§  ,SU(x,T))  +  [VSU]T0  (f  ,SU)  +  (f  ,Sz) 

-  (%w)  +  ((%s  -  S)  v’’sui)  -  ~ 

[(AijViJUj)]  \%z%ldt  +  (L  +  (V\N(U\z')))  \ „TSt  =  0 


(2.68) 


Grouping  the  terms  that  multiply  5U  inside  the  integral  gives  the  adjoint  system  of  PDEs  while 
grouping  the  terms  multiplying  Sz  gives  the  optimal  control  condition.  Notice  the  first  term  in 
the  second  line  multiplies  SUj  not  5Ui  but  since  it  is  summed  over  in  that  term,  any  index  suffice 
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such  that  it  is  legal  to  replace  the  indices  with  j  so  that  all  terms  multiply  the  same  component 
of  SU.  Grouping  the  like  variations  in  SU  at  the  initial  and  final  times  gives  necessary  condition 
on  the  optimal  adjoint  variables. 


J«{§-V',Sz(x,t))  + 


9\  *  .  fdAjj  dUk  _  dAjj  j  y*  _  \  9Vi  _  9Sj  t/*  _j_  dL_ 


’_i _ l  ( 

dt  \  dUk  dx 


dx 


l3  dx  dUi 


V*  +  SU-  + 

vi  '  dUn  ’  UU 3  ' 


(2.69) 


-i  [(2tyV,«yi  ISK  +  (i  +  (V-,N(U\Z-)))  t~Tst  =  o 


^  o  r  \ 

—  -y*(x,0),(J£/(a;,0)j  =0 


(2.70) 


(^  +  v*(x,T),5U(x,T)j  =0  (2.71) 

Lastly,  the  final  time  may  be  free  to  vary.  If  so,  a  necessary  condition  for  the  final  time  must  be 
constructed,  known  as  the  transversality  condition.  The  variation  of  the  state  at  the  final  time 
will  have  a  first  order  term  expanded  in  space  SU (x,  T )  =  SU  +  mt^T>  St.  Therefore,  there  is 
a  change  to  the  final  time  boundary  condition. 


(sU(x,T)  + 


dU(x,  T) 
dt 


St, 


d  K 

dll 


(2.72) 


The  first  term  recovers  the  same  result  as  for  the  fixed  final  time.  The  second  term  remains 
and  is  added  to  another  contribution  to  the  variation  of  the  final  time  from  the  integrand  of  the 
augmented  cost  functional.  This  term  arises  as  a  result  of  Leibniz’s  rule  for  a  derivative  of  an 
integral  with  a  moving  boundary.  When  the  final  time  was  fixed,  the  derivative  with  respect 
to  e  could  be  brought  inside  the  integral  with  no  other  consideration.  Since  the  final  time  (the 
boundary  of  the  time  integral),  is  free  to  vary,  Leibniz’s  rule  dictates  that  the  derivative  can  be 
brought  inside  the  integral  with  the  addition  of  boundary  terms.  For  a  general  function  of  time 
/  over  the  time  dependent  interval  | a (7),  6(f)]  this  looks  like: 
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b{t)  df 

-y,  '  dt  +  f(b(t))b  -  f{a(t))a 
it)  dt 


(2.73) 


d 

dt 


rb(t) 
I  a(t) 


f(t)dt  = 


In  this  problem,  only  the  final  time  is  free  so  a  =  0  and  the  analogous  term  at  bit)  is, 


(L  +  {V,  N(U,  z)))  \t=Tdt  (2.74) 

Combining  both  terms  that  multiply  St  gives  the  inner-product: 

(*• (i  +  l‘-T  +  ff  +  (V, f ))  =  0  (2.75) 

Recall  that  N  =  0  meaning  that  term  contributes  nothing  to  the  sum  at  any  moment  in  time. 
The  term  involving  K  can  be  written  as  a  total  derivative  in  time. 


=  0 


The  Hamiltonian  is  defined  in  general  in  Equation  2.77. 


(2.76) 


H 


L  + 


(2.77) 


Using  this  definition,  the  Transversality  Condition  takes  it’s  familiar  form  in  Equation  2.78. 


+  =0  (2.78) 

Equations  (2.69)-(2.71)  and  (2.78)  gives  a  set  of  first-order  necessary  conditions  that  an  op¬ 
timal  system  (z*,  U* ,V* ,T*)  must  satisfy  if  it  can  be  a  candidate  for  a  local  minimizer  of  J 
constrained  to  N(U,  z ),  the  general  PDE  system  under  a  control  action.  Initial  and  boundary 
conditions  are  either  specified  or  penalized.  Formal  sufficient  conditions  for  optimality  involve 
going  to  second-order  however  it  will  be  obvious  by  the  sequence  of  iterations  of  a  solution 
procedure  if  J  is  being  locally  minimized  or  not. 
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2.3.2  Formulation  of  Single-Phase  Optimal  Control  Shock  Attenuation 
Problem 

This  single-phase  control  calculation  is  meant  to  give  insight  into  a  multi-phase  control  calcu¬ 
lation  where  the  initial  water  placement  determines  the  degree  of  shock  attenuation.  As  stated 
in  the  Introduction,  several  interaction  mechanisms  between  the  water  droplets  and  the  IOP 
wave  are  present  and  not  fully  understood.  The  dominant  dissipative  mechanism  for  the  shocks 
of  interest  is  the  loss  of  energy  of  the  gas  through  vaporization  of  the  water  droplets  [7].  Ex¬ 
perimental  data  from  droplet-shock  interaction  in  this  regime  shows  that  the  other  dissipative 
mechanisms  e.g.  drag  on  droplets,  sensible  heating  of  droplets  etc.  are  less  significant  to  IOP 
attenuation.  In  a  multi-phase  calculation,  the  control  action  would  take  the  form  of  a  liquid 
mass  source  and  the  effect  of  vaporization  will  take  energy  out  of  the  gas  phase.  To  most  sim¬ 
ply  replicate  the  dominant  dissipative  mechanism  with  a  single-phase  calculation,  the  control 
will  act  as  an  energy  sink  distributed  in  space  and  time  behind  the  moving  shock  front.  No 
sinks  or  sources  in  the  momentum  or  mass  conservation  equations  is  needed  for  this  simplified 
formulation. 

The  cost  functional  in  this  context  must  reflect  that  a  decrease  in  the  maximum  jump  in  pressure 
(the  overpressure  at  the  shock  front)  is  most  desirable.  It  should  also  penalize  control  action  but 
to  a  lesser  magnitude.  Therefore,  let  J,  the  cost  functional  be  defined  in  Equation  2.79. 

J——[  [  z(x,  t)2  dxdt  -\ —  /  (P(x,  T)  —  Q(x))2+dx  (2.79) 

2  Jo  Jn  2  Jna 


The  final  time  T,  is  not  fixed,  z(x,  t )  is  the  control  action,  P(x,  T )  is  the  Pressure  at  the  final 
time,  Q(x)  is  the  desired  final  pressure  and  a  and  b  are  weighting  constants.  The  larger  b 
is  compared  to  a,  the  more  significant  the  final  time  penalty  and  the  less  the  penalty  for  using 
control  action.  From  the  form  of  Equation  2.56,  the  individual  cost  functionals  can  be  discerned. 
There  is  no  free  initial  data  and  therefore  no  need  to  associate  a  cost  to  it,  making  1  =  0  for 
the  single-phase  formulation.  The  running  cost  L(U,  z)  need  only  penalize  control  effort  since 
the  final  time  target  state  sufficiently  describes  the  degree  to  which  a  controlled  solution  had  a 
desirable  outcome. 


31 


L(U,z)— —  f  z(x,t)2dx 
2  Jn 


(2.80) 


K(U(;T))J~!  (P(x,T)  -  Q(x))\dx 
2  Jna 


(2.81) 


The  control  action  will  take  the  form  of  an  internal  energy  sink  distributed  in  space  and  time 
behind  the  shock  front.  It  only  appears  on  the  right  hand  side  of  the  energy  balance  equation  of 
the  ID  Euler  system. 


U  =  (p(x,  t),pu(x,  t),pE(x,  t))T 


(2.82) 


d 

{  p  1 

d 

(  pu  \ 

f  °  1 

dt 

pu 

dx 

pu 2  +  P 

— 

0 

\PE  j 

\  u  (pE  +  P)  ) 

^  z{x,t)  y 

pE 


pe  + 


pu 2 

~Y 


(2.83) 


(2.84) 


P  =  pe{  7-1) 


(2.85) 


Initial  conditions  are  fixed  as  stationary,  ambient  air.  Stating  the  density,  velocity  and  pressure 
determines  the  internal  and  total  energy  and  the  conservative  vector. 
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p(x,  0) 

=  1  kg/m 3 

u(x, 0) 

=  0  m/s 

P(x,  0) 

105  Pa 

pu(x,  0) 

=  0  kg/m2  ■  s 

pe{x,  0) 

=  250000  J 

pE(x,  0) 

=  250000  J 

Therefore,  the  initial  variation  in  the  state  vector  in  Equation  2.70  is  zero.  In  addition,  there  is 
no  need  for  an  initial  penalty  I  in  this  formulation  so  the  derivative  =  0  and  V (x.  0)  is  left 
undetermined. 

Remark  on  applying  the  inlet  boundary  condition: 

The  inlet  boundary  condition  is  explicitly  given  by  the  IOP  simulation  data  when  the  flow  is 
supersonic.  If  the  flow  is  subsonic,  non-reflection  of  the  u  —  c  characteristic  is  imposed  [12]. 
In  addition,  the  Monitor  Point  data  was  chosen  where  the  flow  was  nearest  to  ID;  however,  the 
2D  data  did  have  transverse  motion.  The  data  will  still  give  a  plausible  ID  blast  wave  with  the 
inlet  boundary  condition  set  in  this  manner  and  the  goal  of  the  calculation,  controlling  a  range 
of  blast  waves,  can  be  achieved.  At  the  outlet  boundary  the  flow  remains  stationary  because 
the  final  time  will  always  be  such  that  the  shock  wave  will  not  have  enough  time  to  propagate 
though  the  entire  domain  and  reach  the  outlet. 

In  order  to  determine  the  optimal  control  z*(x,t )  that  minimizes  the  cost  functional  J  it  is 
necessary  to  define  the  Hamiltonian  of  the  system  and  derive  necessary  conditions  using  the 
Pontryagin  Minimum  Principle  and  the  calculus  of  variations.  Recalling  Equation  2.77,  the 
Hamiltonian  for  this  system  is  given  in  Equation  2.87. 


H(U,V,z,t ) 


f  a  2  dp  d  (pu)  d  (pE) 

In  P  +Vlm+v=—  +  V3—dx 


(2.87) 
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T 

The  co-state  or  adjoint,  vector  is  (Vj,  V2,  V3)  .  The  flow  dynamics  will  be  solved  in  the  con¬ 
servative  form  of  Equation  2.83,  yet  it  is  useful  to  write  the  system  in  the  quasi-linear  form  of 
Equation  2.88  to  derive  the  adjoint  system.  The  conservative  basis  Jacobian  matrix  was  given 
in  Equation  2.7.  Now  the  system  has  a  control  source  vector  in  the  energy  balance  equation. 


dUi 

dt 


+  Aij(U) 


Wi 

dx 


^i^i3 


(2.88) 


To  derive  necessary  conditions  as  in  the  general  case,  the  optimal  state  must  be  defined 
(U*(x,  t),  V*(x,  t),z*(x,  t),T*)  and  the  optimal  control  perturbed  such  that  z  =  z*  +  e8z  where 
e>  0  is  a  small  constant.  The  variation  in  the  control  will  cause  variational  terms  in  each  of 
the  other  free  variables  of  the  system  that  must  necessarily  vanish  at  an  optimal  solution.  To 
incorporate  the  constraints  of  the  ID  Euler  system  using  the  Lagrange  multiplier  method,  each 
conservation  law  is  multiplied  by  an  adjoint  variable  and  added  to  J,  resulting  in  the  augmented 
cost  functional  J.  Then  from  expanding  J  in  a  Taylor  series  the  first  order  necessary  condition 
will  be: 


-7- J  (z*  +  eSz)  |e=0  =  0  (2.89) 

de 

Grouping  the  terms  of  like  variational  multipliers  gives  the  optimal  system.  Integrating  by 
parts  until  all  derivatives  are  on  the  adjoint  vector  yields  the  following  linear  system  of  partial 
differential  equations  [11]. 


^  4.  A.  4.  (™*L 
dt  VJ  dx  \  dx 


dAij 

dt 4 


c)  c)  L 


0 


(2.90) 


Non-trivial  elements  of  the  matrix 


dAij 

9Uk 


are  given  in  Equations  2.91-2.95. 


d  A  dUk 
dUk'  dx 


(7-3) 


u2  dp 
p  dx 


u d ( pu ) 
p  dx 


(2.91) 
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(2.92) 


d  dUk 
dUk'  dx 


(7-3) 


f  u  dp 
\p  dx 


1  9  (pu)  \ 
p  dx  ) 


d  A  dUk 
dUk'  dx 


~p  -  3  (7  -  1)  «2)  ^  ^  (3  (7  -  “2 


d  ( pu ) 
dx 


u  d  ( pE ) 

^  p  dx 

(2.93) 


d  dUk 
dUk'  dx 


i(3(T-l)U2-7s)g 


3(7-1) 


u  d  (pu) 
p  dx 


,  7  8<J>E) 
p  dx 


(2.94) 


_d_A  dUk  _  dp  7  d  (pu) 
dUk  33  dx  p  dx  p  dx 

dQ  denotes  the  boundary  of  the  spatial  domain.  The  right  hand  side  of  Equation  2.90  is  zero  in 
this  formulation  because  the  running  cost  does  not  explicitly  depend  on  the  state  vector. 

Since  the  final  state  is  not  fixed,  there  is  a  necessary  condition  on  the  adjoint  vector  at  the  final 
time. 


V?(z,n=  -^jK(U(x,r)) 


Therefore,  in  this  basis,  all  3  derivatives  are  non-zero  in  f2s. 


(2.96) 
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V\  (x,  T*)  =  b{P{x,T*)-Q(x))+~{x,T*) 

=  b{P{x,T*)-Q{x))+-U{x'P2  1) 


V*(x,T*)  =  b{P(x,T*)-Q(x))+.-^r(x,T*) 

-b  (P  {x,  T*)  -  Q(x))+  ■  u  {x,  T*)  (7  -  1) 


(2.97) 


V3*{x,T*)  =  b  (P  (x,  T*)  —  Q(x))_ 


dP 


d{pE) 

=  b  (P  (x,  T*)  —  Q(x))+  •  (7  —  1) 


(x,T*) 


For  a  free  final  time,  the  necessary  condition  for  the  optimal  final  time  T*  is  given  by  the 
Transversality  condition.  Define  /  as  a  functional  to  be  used  in  the  solution  procedure. 


H  (U*  (x,  T*) ,  V*  (x,  T*) ,  z*  (x,  T*) ,  T*)  +  —K  (U  (x,  T*))  = 

at 


/  %z(x,  T*)2  +  b(P(x,  T*)  -  Q(x)) 
Jn  2 


dx  = 


=  f(T* 


(2.98) 


The  necessary  condition  on  the  optimal  control  solution  comes  from  maximizing  the  Hamilto- 
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nian  for  an  unconstrained  control.  The  unconstrained  condition  is  justified  because  there  is  no 
restriction  on  control  magnitude  in  regions  where  control  is  allowed.  The  integral  is  true  over 
any  domain  Q  and  at  any  moment  in  time  it  should  be  true  point- wise  for  all  t  G  [0,  T], 

r 

——([/*,  V*,  z*,  t)  —  /  az*(x,  t)  +  14* (x,  t)dx  —  0  (2.99) 

oz  Jn 

Equations  (2.82)-(2.86),(2.90),(2.97),(2.98)  and  (2.99)  give  the  complete  set  of  first  order  nec¬ 
essary  conditions  for  the  optimal  system. 

2.3.3  Formulation  of  Two-Phase  Optimal  Control  Shock  Attenuation  Prob¬ 
lem 

In  the  two-phase  system,  the  control  action  takes  the  form  of  the  free  initial  data.  The  first 
state  variable  oti(x,  0),  the  initial  water  volume  fraction  distribution,  will  solely  determine  the 
attenuation  at  the  final  time.  This  is  a  satisfying  formulation  as  the  alternative  would  involve 
water  injection  over  time  as  the  control.  However,  water  cannot  be  significantly  injected  into  a 
flow  over  the  interaction  time  scales  of  a  shock  wave.  The  goal  of  attenuating  the  jump  in  gas 
pressure  at  the  final  time  is  unchanged,  so  the  final  time  penalty  need  only  adopt  the  two-phase 
notation.  Equation  2.100  gives  the  cost  functional  for  the  two-phase  system. 

J=  [  ”  (cn(x,  0))2  dx  +  f  ^  (pg(x,T)  -  Q(x))2+dx  (2.100) 

Note  that  there  is  no  running  penalty,  thus  L  =  0  in  this  formulation.  It  is  then  clear  what  the 
initial  penalty  and  final  time  penalties  are,  shown  in  Equations  2.101  and  102  respectively. 

I(U(;  0))  =  (on(x,  0))2  dx  (2.101) 

K(U(;T))=  f  *(Pg(x,T)-Q(x))2+dx  (2.102) 

The  multiphase  system  can  be  defined  in  the  primitive  basis  recalled  from  Equation  2.32. 
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(2.103) 


U  (dig,  Pg,  Ug,  Pg,  Pl,  Ul,  P[) 

The  state  dynamics  are  of  the  form  of  Equation  2.34  with  the  Jacobian  matrix  defined  in  Equa¬ 
tion  2.33.  The  adjoint  system  has  the  form  of  Equation  2.69.  Note  the  inclusion  of  the  term 

present  in  the  two-phase  system  but  not  the  single-phase  system.  This  term  as  well  as  other 
non-trivial  closed  form  elements  of  the  matrices  of  the  adjoint  system  are  given  in  Appendix  A 
and  B  for  brevity.  In  the  primitive  basis,  the  density,  velocity  and  pressure  of  the  gas  and  liquid 
have  fixed  initial  conditions  which  are  stationary  ambient  atmosphere  at  sea  level. 

pg(x,  0)  =  1  kg/m 3  pi(x,  0)  =  103  kg/m 3 

ug(x,  0)  =  0  m/s  ui(x,  0)  =  0  m/s  (2.104) 

Pff(x,0)  =  105  Pa  Pt(x,  0)  =  105  Pa 

The  first  variable,  the  liquid  volume  fraction,  is  not  fixed  at  the  initial  time  and  acts  as  the 
control.  The  initial  cost  functional  /  is  then  associated  with  a  penalty  to  the  amount  of  water 
used. 

Looking  back  to  Equation  2.70  and  combing  Equation  2.101  and  2.105  gives  Equation  2.106. 

<J£7i(a;,0)  =  /  °  (2.105) 

(^’5^(X’°))  =  (2.106) 

After  dropping  the  factor  of  6ag(x ,  0)  the  first  order  necessary  optimal  condition  on  V/(x.  0)  is 
given  in  Equation  106. 

[  a  (1  —  ag(x,  0))  —  V*(x,  0)dx  —  0  (2.107) 

Since  this  is  true  over  a  general  domain  f)  it  must  be  true  point-wise. 

At  the  final  time  T,  our  cost  functional  K  will  penalize  the  gas  pressure  profile  in  space  if  the 
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overpressure  at  the  shock  front  is  greater  than  that  of  a  predefined  target  state.  K  is  only  an 
explicit  function  of  the  fourth  state  variable  f/4,  the  gas  pressure.  The  rest  of  the  state  variables 
are  free  to  vary  and  their  variation  does  not  vanish  at  the  final  time. 

{0  i  +  4 

,  .  .  (2.108) 

7^0  i  =  4 

Using  the  condition  in  Equation  2.71  with  the  definitions  in  Equations  2.102,  2.108  and  2.109 

gives  final  time  necessary  condition  in  Equation  2.110. 

(H, «/<(*, T))  =  SP^x.T))  (2.109) 

Therefore,  as  before,  the  final  time  optimal  condition  on  V£(x,  T )  becomes  Equation  2.110: 

(H :-  +  v:(x,T),8U4(x,T)j  =0  (2.110) 

where  Q(x)  is  the  predefined  target  state  and  b  is  a  constant.  The  optimal  condition  then  be¬ 
comes  Equation  2.111: 


I  b{Pg{x,T)  ~  Q(x))  +  V;(x,T)dx  =  0  (2.111) 

J  Q 

which  must  also  be  true  point-wise. 

The  purpose  of  the  calculation  is  to  optimize  blast  wave  attenuation.  To  that  end,  let  the  forcing 
function  or  blast  wave  condition  at  a  stationary  point  over  time  be  the  left  boundary  condi¬ 
tion  UL(t)  with  a  right- ward  traveling  blast  wave.  The  flow  conditions  are  supersonic  and  may 
be  considered  fixed  unless  reflected  shock  waves  originating  from  within  the  simulation  do¬ 
main  must  propagate  outward  in  a  non-reflection  characteristic  boundary  condition.  The  spatial 
boundary  term  from  Equation  2.69  for  a  linear  interval  x  €  [xl,xr]  =  Q  is  simplified  in 
Equation  2.1 12. 
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£  IK*’ 5(l)]  =  £jn  AV*  ■  SUdx  =  /:=:«  £  (AV*  •  SU)  dx  =  [AV*  ■  SU]^* 

=  A  ( UR )  VS  •  -  A  (t/L)  VL  •  Sill 

=  A(u (xr,  t))  V*(a;fl,  t)  •  5*7 (xfl,  t)  —  A(u (xl,  t))  V*(xL,  t)  •  SU (xL,  t ) 

(2.112) 

The  calculation  can  be  performed  in  such  a  way  that,  after  starting  with  no  motion  in  the  initial 
conditions,  the  blast  wave  can  be  simulated  as  entering  at  the  left  boundary,  moving  right- ward, 
and  not  reaching  the  right  boundary  by  the  final  time.  In  such  a  simulation  there  is  no  change 
to  the  state  vector  at  the  right  boundary  for  all  time. 

{0  X  =  Xr 

0  x  =  xl  supersonic  inflow  (2.113) 

7^  0  x  =  xl  subsonic  inflow 

From  Equation  2.112,  the  boundary  integral  from  Equation  2.111  reduces  to: 


40 


Because  we  have  all  the  control  effort  taking  place  at  the  initial  time,  through  the  initial  data, 
there  is  no  need  for  control  over  time.  So,  z*(x,t )  =  0  and  since  it  is  fixed,  the  variation 
5zi(x,  t)  vanishes  as  well. 

The  final  time  penalty  K  can  completely  communicate  the  objective  of  overpressure  attenuation 
without  the  need  for  a  penalty  over  space  and  time  L  for  matching  a  target  over  space  and  time. 
All  of  this  means  that  L(U,  z)  —  0  and  consequently,  the  first  two  inner  products  inside  the  time 
integral  in  Equation  2.68  are  also  zero.  The  complete  set  of  first  order  necessary  conditions 
for  the  two-phase  optimal  system  is  given  by  Equations  2.103-2.104,  2.107,  2.111,  2.115  given 
the  state  dynamics  of  Equation  2.34  and  corresponding  adjoint  PDE  from  Equation  2.69.  The 
closed  form  matrix  elements  for  the  adjoint  PDE  are  stated  in  Appendices  A  and  B. 
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CHAPTER  3: 
Numerical  Methods 


3.1  Computational  Fluid  Dynamics  for  Unsteady  Compress¬ 
ible  Flow 

3.1.1  Exact  Solution  to  the  Riemann  Problem 

This  section  describes  a  numerical  method  [40],  [57]  used  for  calculating  the  exact  solution 
to  the  Riemann  Problem  (Equation  2.16)  for  the  Unsteady  Euler  Equations  of  gas  dynamics 
(Equations  2. 1-2.5).  The  method  takes  advantage  of  the  knowledge  of  the  wave  patterns  shown 
in  Figures  2.1  and  2.2.  The  subscripts  L  and  R  refer  to  the  left  and  right  initial  data  and  constant 
states  to  the  left  UL  =  ( Pl,ul ,  Pl )  and  right  UR  =  (pR,  uR,  pR)  of  the  star  region  respectively. 
Originating  from  x  =  0  at  an  initial  moment  in  time,  a  rarefaction  wave  moves  to  the  left  at 
speed  u  —  c,  a  contact  discontinuity  moves  to  the  right  at  a  speed  u  and  a  shock  wave  moves  to 
the  right  at  a  speed  u  +  c.  The  speed  of  sound  changes  across  the  star  region  making  it  locally 
determined.  The  exact  solution  is  constructed  by  sampling  the  flow  at  any  moment  in  time  t  >  0 
on  a  discrete  grid  and  iteratively  solving  Equation  3.1  using  knowledge  of  the  sample  points’ 
location  in  relation  to  the  wave  pattern  and  given  the  definitions  in  Equations  3. 2-3. 4  previously 
stated  in  Chapter  2. 

/  (p\  UL,  UR)  =  fL  (p,  Ul)  +  In  (p,  Ur)  +UR-Ul  (3.1) 


(3.2) 


(3.3) 
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2 


AL  = 
Ar  = 
Bl  = 
Br  = 


(7+1  )Pl 


(7+1  )PR 

U~ l\T 

(7+l)PL 
(7-1) 


(7+1) 


PR 


(3.4) 


For  each  discrete  sample  point,  the  process  starts  with  an  initial  guess  for  the  pressure  p*.  Next, 
the  value  of  f(p*,  Ul,  Ur)  is  determined  based  on  the  known  wave  patterns.  The  velocity  u*  is 
updated  according  to  Equation  3.5  or  3.6  depending  on  the  part  of  the  star  region  the  sampling 
point  is  located. 


u*L  =  uL-  fL(p*,UL)  (3.5) 

u*r  =  uR  +  fR(p*,UR)  (3.6) 

Next,  update  p*  according  to  Equation  3.7  or  3.8. 

=  1/7  (3-7) 

P‘R  =  PR  (^f)  <3'8) 

The  shock  speed  in  a  stationary  frame  is  denoted  as  S.  This  solution  method  is  impractical  since 
many  iterations  are  required  to  solve  for  the  root  of  Equation  3. 1  for  each  discrete  point  in  space 
at  each  instant  in  time.  The  exact  solution  is  shown  in  red  in  Figure  3.1. 

3.1.2  Godunov  Method  with  an  Approximate  Riemann  Solver 

The  first-order  Godunov  Method  [58]  assumes  a  piece-wise  constant  solution  in  each  discrete 
cell  volume.  The  Riemann  Problem  is  then  solved  locally  between  each  cell  at  each  moment 
in  time  based  on  the  previous  piece-wise  solution.  This  approximate  scheme  does  not  require 
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an  iterative  solution  and  therefore  is  much  less  computationally  demanding.  With  the  choice 
of  a  conservative  method,  the  errors  introduced  by  the  approximate  solution  behave  in  a  known 
and  acceptable  way  [59]-[64].  Figure  3.1  compares  approximate  numerical  solutions  of  the 
shock  tube  problem  to  the  exact  solution,  four-tenths  of  a  millisecond  in  time  after  the  di¬ 
aphragm  bursts.  The  finite  difference,  Modified  MacCormick  method  [65]  introduces  dispersive 
error  while  the  conservative  Godunov  method  introduces  dissipative  error.  The  dispersive  error 
causes  the  wrong  jump  condition  and  therefore  the  wrong  entropy  solution,  shock  speed  and 
position  as  shown  in  Figure  3.1.  Similar  results  exist  in  the  literature  [66].  More  oscillations 
in  the  solution  will  also  cause  more  problems  when  equilibrating  with  another  phase.  On  the 
other  hand,  dissipative  errors  are  second  order  in  spatial  derivatives  of  U  and  therefore,  are  a 
viscous-like  error.  It  is  well  established  that  in  the  limit  of  vanishing  viscosity  the  correct  shock 
solution  is  achieved.  From  a  physical  point  of  view,  a  truly  inviscid  flow  is  being  modeled. 


xIO5 


Figure  3.1:  Comparison  of  numerical  solutions  to  the  Riemann  problem  .4  ms  after  diaphragm 
burst  in  shock  tube  flow.  The  initial  pressure  ratio  was  10:1  in  each  case  with  a  discontinuity  at 

x  =  .5  m 


Let  the  spatial  domain  x  G  [xL,xR]  =  Q  be  divided  into  M  intervals  each  equally  spaced  by 
Ax  =  ( xR  —  xL)  /M.  The  locations  of  the  discrete  grid  points  are  given  in  Equation  3.9,  for 
j=l,...,M+l. 
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Xj  =  (j  —  l)Aa; 


(3.9) 


The  solution  will  be  piece-wise  constant  between  half-cells.  As  a  sensible  choice  for  limits  of 
integration,  the  half-cell  locations  are  define  in  Equation  3.10,  for  j  =  1,...,M. 


Xj- 1/2  =  (j  -  1/2) Ax  (3.10) 

For  all  cell  volumes  let  the  limits  of  integration  be  x\  =  Xj-1/2  and  X2  =  Xj+ 1/2.  Let  the  time 
domain  [0,T]  be  divided  into  N  intervals  each  with  size  At  =  T /N.  The  moments  in  time  are 
given  in  Equation  3.11  for  n  =  1,...,N. 


tn  =  {n-l)At  (3.11) 

For  all  regions  of  the  discrete  space-time  domain  [xj-i/2,  ^+1/2]  x  [tn,tn+1]  the  conservation 
law  in  Equation  2.21  becomes  Equation  3.13  without  approximation. 


f  F(U (xj-1/2,  t))dt  —  j  F(U(xj+1/2lt))dt 

J  t\  J  t\ 

(3.12) 


Discrete  approximations  to  the  solution  at  two  discrete  moments  in  time  are  defined  in  Equation 
3.13  and  3.14. 


1  — 

U (x,  tn)dx  (3.13) 

Ax 

1  i-\-l /2  — 

U?+1^—  U (x,  tn+1)dx  (3.14) 

Ax  Jxj_1/2 

The  cell  average  of  the  conservative  vector  is  denoted  as  U(x,t).  Locally  in  each  cell,  U  is 
constant  when  traveling  at  the  estimated  wave  speeds.  This  is  a  valid  assumption  as  long  as 
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the  time  step  is  sufficiently  small  such  that  no  wave  is  propagated  beyond  the  cell  from  where 
the  speed  was  estimated  to  avoid  wave  interaction  which  causes  inaccuracy.  This  leads  to  the 
famous  CFL  condition  given  in  Equation  3.15. 


At  < 


Ax 

IT 


(3.15) 


This  time  interval  is  limited  by  the  maximum  wave  speed  in  each  cell  S.  For  t  G  [tn,  tn  +  At] 
the  discrete  inter-cell  solution  and  flux  are  given  in  Equations  3.16  and  3.17. 


Uj-i/2  —  U(xj^il2,tn)  (3.16) 

Fj_1/2  =  F{u(xj^1/2,tn))  (3-17) 

With  the  discretizations  of  Equation  3.13,  3.14,  3.16  and  3.17  and  the  CFL  condition  on  At,  the 
integral  form  of  the  conservation  law  in  Equation  3.12  gives  the  numeric  form  of  Godunov’s 
Method  in  Equation  3.18. 


Ur'  =  U?-A(Fj+1/2-Fj-1/2)  (3.18) 

A  suitable,  conservative  approximation  to  the  inter-cell  flux  FJ+1/2  will  result  exclusively  in 
dissipative,  not  dispersive,  numerical  error  near  the  shock  front.  As  described  earlier,  these 
attributes  are  necessary  for  accurate  jump  conditions  at  the  shock.  In  the  limit  of  vanishing 
viscosity  they  converge  to  the  correct  weak  solution.  The  Rusanov  approximation  used  for  the 
inter-cell  flux  [40],  [54]  is  given  in  Equation  3.19. 

Fj+ 1/2  =  \(F  ( U y  +  F  (Uj+1))  -  S,  (Ui+l  -  Uj)  (3.19) 

The  maximum  wave  speed  in  the  jth  cell  is  S3 .  For  a  shock  this  is  given  by  the  Rankine- 
Hugoniot  jump  condition  in  Equation  2.13  of  Chapter  2. 


47 


3.1.3  Second  Order  Godunov  Method 


A  conservative  second  order  extension  to  the  Godunov  method  [67]-[70]  was  implemented  to 
solve  the  compressible  flow  dynamics  in  ID  under  a  distributed  control  action.  This  method  as¬ 
sumes  that  the  solution  in  each  cell  is  piece- wise  linear  and  projects  an  intermediate  solution  on 
a  non-uniform  grid  based  on  the  maximum  characteristic  speeds  from  the  interpolated  solutions 
at  each  cell  interface.  Figure  3.2  is  a  diagram  of  what  the  uniform  and  non-uniform  grids  might 
look  like,  based  on  the  wave  speed  estimates  at  the  inter-cell  boundary. 


Figure  3.2:  Intermediate  non-uniform  grid,  Kurgarov  and  Tadmor  [67] 


The  familiar  Godunov  integration  on  the  uniform  grid  is  then  second  order  accurate  because  of 
the  intermediate  finer  grid.  This  numeric  method  was  chosen  because  it  preserves  the  conserva¬ 
tive  property  with  only  dissipative  error  and  can  be  adapted  to  the  solution  of  a  ID  multiphase 
calculation  as  presented  in  Saurel  and  Abgrall  [31].  The  essential  numerical  steps  are  summa- 
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rized  here. 


At  a  current  time  step  tn  the  solution  vector  [/"  is  known  for  all  x3  .  The  local  spatial  derivative 
of  the  solution  vector  (Ux)n  can  be  calculated  by: 


minmod  (  3 


U?- 1 


TTTl 

Uj  + 1 


Ax 


Ax 


(3.20) 


The  minmod  function  provides  an  up  wind/ down  wind  switch  based  on  the  wave  direction  of 
motion. 


minmod  (a,  b ) 


-  ( sign  (a)  +  sign  ( b ))  •  min  (|a|,  |6|) 


(3.21) 


The  strategy  is  to  first  create  a  finer,  intermediate  grid  to  calculate  the  solution  iv”  1 1/2, 
at  a  half  time  step  forward  in  time  with  a  conservative  Godunov  method.  The  second  stage  is 
to  calculate  the  solution  U!-'+ 1  at  the  next  time  step  on  the  coarser,  uniform  grid  with  a  conser¬ 
vative  Godunov  method.  The  result  is  that  spatial  resolution  at  the  shock  front  more  accurately 
captures  the  shock  front  discontinuity  while  maintaining  a  conservative  second  order  dissipa¬ 
tive  error  shape.  Other  methods  can  introduce  dispersion,  losing  the  conservative  property  and 
consequently  will  have  oscillations  near  the  shock  front  which  give  the  wrong  jump  conditions. 

The  upper  bound  on  the  maximum  wave  speed  at  the  inter  cell  boundary  as  estimated  with  a 
first  order  Taylor  expansion  of  w” . 


«"+ 1/2  =  ™x  ( eig  (A  (Uj+1/2))  ,  eig  [A  (u++1/2)))  (3.22) 


77+  —  TTn 

Uj+ 1/2  -  Uj+ 1 


Ax 


(Us 


n 

xJj+ 1 


(3.23) 


U~+l/2  =  f/’Vi  +  ~y  (tf*)"  (3.24-) 

The  superscripts  +  and  ~  refer  to  the  solution  to  the  right  and  left  of  the  inter-cell  discon¬ 
tinuity.  The  eigenvalues  of  A  are  the  eigenvalues  of  the  characteristic  wave  speed  es- 
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timated  at  the  inter-cell  boundaries.  For  the  single-phase  Euler  system  the  eigenvalues  are 

(ug,  ug  +  Cg,  Ug  —  Cg )  and  for  the  two-phase  system  are  (V*,  ug ,  ug  —  cg,  ug  +  cg,  ui,  ui  —  ci,  ui  +  q). 

The  finer,  non-uniform  intermediate  grid  is  defined  by  using  the  bound  on  the  maximum  wave 
speed  at  each  inter  cell  boundary. 


—  rrn  —  nn  ■  At 

xj+l/2,l  ~  xj+ 1/2  ajf+l/2  L*-i 


(3.25) 


Xj+l/2,r  —  x^+ 1/2  +  aj+l/2  '  At 


(3.26) 


The  solution  on  the  non-uniform  grid  is  given  in  Equations  3.27  and  3.28. 


w 


3+ 1/2 


VJ ■  +  I£l  +  A»-.Jfl/aAt( 


9nn 

Zllj+ 1/2 


f  (q+X)  -  f  (x/i) 


n+l/2  \ 


(3.27) 


wj  —  Uj1  +  2  (ai-i/2  1/2)  (^®) 

At/ Ax 


1— At/ Ax 


aj-l/2+aj+l/2 


P  (  rrn+1/2 


P  ( rrn+1/2 
^  X/+l/2,r 


(3.28) 


Quantities  on  the  uniform  grid  at  the  inter-cell  boundaries  are  approximated  at  the  current  time 
step  in  Equation  3.29  and  3.30  and  at  the  half-time  step  in  Equations  3.31  and  3.32. 


q+iW  =  U]  +  A*(ty”  (i  - 


n 

7+1/2 


(3.29) 


q+1/2,.  =  c?+1  +  Ar  (cy;+1  (1  -  |^"+1/2 


(3.30) 


rrn+l/2  _  jjn  _  rr  / TTTI 

Uj+l/2,l  ~  Uj+l/2,l  9  r  (Uj+l/2,l 


(3.31) 
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T  7- 77. +  1/2 
Uj+l/2,r 


U?+  l,2,r~YF(U 


3+1/2, r 


(3.32) 


The  solution  is  progressed  to  tn+l  on  uniform  grid  with  a  conservative  second  order  update  as 
shown  in  Equation  3.33  and  using  the  definition  in  Equation  3.34. 


rrn+l  _  At  n  mn+1  -+  ( 1  —  Af  ( nn  i  nn  ))  i  ALnn  mn+l  4- 

Uj  ~  Axaj-l/2Wj-l/2  +  Ax  \aj-l/2  +  aj+l/2j )  Wj  +  Axaj+l/2W j+l/2^ 


Ax 

2 


(( 


3~  1/2 


n+ 1 
j-1/2 


n+1  3 


(3.33) 


(U: 


\n+l 


x>j+l/2  ~  ^ 


=  - — minmod 


w 

1  +  M( 

1  '  Ax  ( 


n+1 

3+i 


w 


n+1 

3+1/2 


3+1/2  Uj+ 3/2 


W 

At 


n+1 

3+1/2 


w 


n+1 


i +  —  r 

1  ^  As  V 


3+1/2  Uj+3/2j 


(3.34) 


3.2  Method  of  Solution  for  Two-Phase,  Gas-Liquid  Coupling 

The  two-phase  balance  equations  from  Equations  2.22-2.24,  the  point-wise  relations  given  in 
Equations  2.25-2.28,  the  interface  quantities  defined  in  Equations  2.29-2.31  and  source  terms 
defined  in  Equations  2.35  of  Chapter  2  are  the  basis  for  the  overarching  numeric  method  used 
to  couple  the  two-phase  system  based  on  the  model  from  Saurel  and  Ab grail  [31].  A  survey  of 
the  literature  on  numeric  methods  for  two-phase  flow  [71]-[76]  reinforces  the  advantages  of  this 
method  and  the  theory  on  which  it  is  based. 

3.2.1  Hyperbolic  Operator 

As  was  stated  in  Chapter  2,  both  phases  are  compressible,  which  produces  seven  distinct  wave 
speeds  in  the  system.  Equations  2.23-2.24  can  then  be  factorized  into  a  hyperbolic  PDE  with 
a  source  term  H (U)  from  the  motion  of  the  interface  and  three  separate  ODEs  for  the  separate 
components  of  the  non-conservative  source  term  S(U)  [31]. 

+  035) 
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(3.36) 


S(U) 


Ls  =  LMh  •  LVr  ■  Lpr  (3.37) 

When  the  first-order  Godunov  scheme  is  used,  factorization  of  the  two-phase  model  occurs 
according  to  Equation  3.39. 

Un+i  =  LftL^Un  (3.38) 

When  the  second-order  Godunov  scheme  is  used,  factorization  of  the  two-phase  model  occurs 
according  to  Equation  3.40  [76]. 


jjn+l  =  LAt/2LAtLAt/2jjn  (3.39) 

The  ability  of  an  approximate  numeric  method  to  capture  the  shock  discontinuity  is  based  on 
the  approximate  Riemann  fluxes  at  the  inter-cell  boundary.  The  method  gives  a  few  choices  for 
conservative  approximate  Riemann  fluxes.  Equation  3.41  is  the  analogous  Godunov  update  in 
Equation  3.18  for  the  two-phase  system  with  a  first-order  Rusanov  flux  used  to  update  the  mass, 
momentum  and  energy  of  the  gas  and  liquid  phases. 


jjn+l 


Ui- 


A  t 
Ax 


+  AtH  ( U\ *) 


d(Xgi 

dx 


(3.40) 


Fi+ 1/2  —  F  (Ui,  Ui+l) 


^(Fi  +  Fi+1-5'(£7i+1-£7i)) 


(3.41) 


_  ^g,i- (-1  1 

dx  2  ■  Ax 


(3.42) 


an+1  =  an  ■ 

g.i 


At 


2  •  Aa;  L 


9,i+ 1 


—  a 


g,i 


+  S;- 


1/2 


{<i  - 


a 


9,i-l 


(3.43) 
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The  Godunov  method  for  the  two-phase  system  in  Equation  3.41  can  be  extended  to  second 
order  in  an  analogous  was  as  that  described  in  Section  3.1.3. 

3.2.2  Source  Term  Operators 

Next  comes  handling  of  the  source  terms.  They  can  be  factored  into  three  separate  terms  based 
on  timescale  or  phenomenon  and  solved  sequentially  as  three  separate  ODEs  after  the  hyper¬ 
bolic  Godunov  solution.  First  the  pressure  relaxation  operator  in  Equation  3.38  is  isolated  and 
the  solution  is  given  in  Equation  3.44  where  Un  is  the  solution  vector  after  the  Godunov  inte¬ 
gration. 


UPR  =  Un  +  LPR(Un )  (3.44) 

This  ODE  simplifies  to  Equation  2.41  and  2.42  which  are  recalled  below. 

raPR 

(agpgEg)PR  -  ( agpgEg)°  =  /  9  Pidag  (3.45) 

Jao 

{alPiEi)PR  -  (cqp^)°  =  -  9  Pidag  (3.46) 

Ja° 

Internal  energy  is  adjusting  in  each  phase  due  to  the  work  done  by  inter  facial  pressure.  The 
goal  is  to  calculate  aPR  such  that  the  pressures  are  equilibrated  Pg  =  Pi  for  all  Xi  while  also 
satisfying  the  integral  constraints  in  Equations  3.45  and  3.46.  A  solution  procedure  used  in  [31] 
is  summarized  here. 

1.  Make  an  initial  guess  for  aPR  based  on  the  difference  Pg  —  /) 

2.  Compute  ( ag ,  pg ,  Eg)PR  and  («z,  pz,  E{)PR  from  the  integrals  in  Equations  3.46  and  3.47 

3.  During  relaxation,  total  mass  (aff,  pg)  is  constant  so  a  new  iterate  of  pg  can  be  computed 
and  then  Eg  and  eg  as  well.  The  procedure  for  liquid  is  identical. 

4.  Calculate  Pg  and  Pi  from  the  equations  of  state 

5.  Repeat  until  equilibrium  condition  Pg  =  Pi  is  reached  for  all  Xi 
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The  velocity  relaxation,  or  drag,  operator  taking  the  form  of  Equation  3.47  is  next  solved  nu¬ 
merically. 


UVR  =  UPR  +  Lvr(Upr)  (3.47) 

The  drag  force  is  based  on  the  empirical  drag  law  for  spherical  objects  as  given  in  Equation 
2.43  and  shown  in  discrete  form  in  Equation  3.48. 


Fd(xi,  tn)  =  Cdpg(xi,tn)Dl(xi,tn)ai(xi,tn)  ( ug(xhtn )  -ui(xi,tn))2  (3.48) 

With  the  drag  force  on  the  droplet  known,  the  solution  to  the  ODE  is  straightforward,  as  shown 
in  Equation  3.49. 


Ug(Xi,tn+1)  =  Ug{Xi,tn)  + 
Ui(Xi,tn+1)  =  Ui(Xi,tn )  - 


At 


Ctg{Xi,tn)pg(Xi,tn ) 

At 


Fd(xi,tn) 


Eg(xi,tn+1 )  =  Eg(xi,tn )  + 

Ei(xi,tn+1)  =  El(xi,tn)  ~ 


ai(xhtn)pi(xi,tr 

At 


ag(xhtn)pg(xi,tn ) 

At 


Fd[Xi,  tn) 

Fd(xi,  tn)Vi(xi,  tr 


ai(xi,tn)pi(xi,tr 


Fd(xi,tn)Vi(xittr 


eg(xi,tn+1 )  =  eg(xi,tn)  +  -Ug(xi,tn+1Y 
ei(xi,tn+1 )  =  ei(xi,tn )  -  ^ui(xi,tn+1f 


(3.49) 


The  third  and  final  source  vector  left  is  due  to  the  mass  and  heat  exchange  of  the  phases  due  to 
vaporization  and  interface  heat  transfer. 


Un+ 1 


umh  =  uvr  +  Lmh(juvr) 


(3.50) 
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Vaporization  dynamics  of  the  droplets  is  given  by  the  Empirical-Beta  Vaporization  Law  with 
experimental  constants  given  [7].  The  present  calculation  assumes  that  droplets  are  mono- 
dispersed  and  spherical.  Their  initial  size  is  denoted  D®  (xi,  0).  When  surrounded  by  hot  gas, 
the  droplets  will  start  to  vaporize  and  their  diameter  will  evolve  according  to  Equations  3.51 
and  3.52. 


d  D2 

-^(xi,tn+1)  =  -P(Tg(xi,tn)) 


(3.51) 


P  ( Tg(xi ,  tn ))  =  7600  (l  +  7.4  •  10"7  (Tg(xh  tn )  -  30077) +)2'7548  ^m2 /s 


(3.52) 


Once  the  new  size  of  the  droplets  is  known,  the  added  amount  of  gas  volume  fraction  is  known. 


A  ( Dn+l  Dn 3 

an+1  =  an  +  7r  — - T  |  tv" 

»  5  3  2  2  I  p 


(3.53) 


The  critical  quantity  for  sinking  energy  from  the  gas  via  vaporization  is  m  which  is  given  in 
discrete  form  below  based  on  Equation  2.44. 


rh(xi,  tn+1)  =  pH2Ov(Pg(xi ,  tn+1)/ atm)  (ag(xh  tn+1)  -  ag(xh  fn))  (3.54) 

The  density  of  water  vapor  changes  an  order  of  magnitude  in  the  pressure  range  of  1-10  atm.  It 
is  calculated  based  on  a  quadratic  fit  to  steam  table  data  [42]. 


Ph2Ov  ( Pg/ atm )  (Xi  ,tn+1) 


-.0024344- 


p9(xhtn+1y 


atm 


+.5368 


Pg(Xj,tn+1) 

atm 


-.077246  kg/m 3 


(3.55) 

With  the  additional  mass  of  water  vapor  in  the  gas  phase,  the  gas  density  will  increase  according 
to  Equation  3.56  while  the  density  of  liquid  water  does  not  change. 


p+  =  i « 


(«>S  +*)  /« 


?l+l 

9 


(3.56) 
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The  time  integration  can  now  be  performed,  completing  the  numerical  solution  of  the  source 
vector  S(U)  at  each  time  step. 

(agPg)MH  =  Ug+1Png+1  +  ™ 

(agpgug)MH  =  a^+1p'^+1ulgR  +  rhVi 

(agPgEg)MH  =  ang+1png+1EjR  +  m  (. Lhv  +  Et)  +  AtQi 

(  \MH  n+l  n+ 1  •  (3.57) 

(aipt)  =  oil  Pi  ~m 

(■ cqpiUi)MH  =  a?+1pTuYR  -  linVi 

{alPlE{)MH  =  o$+1p?EYR  -  m  (Lhv  +  E>)  -  AtQi 

3.3  Solution  Procedure  for  Optimal  Control  System 

The  method  for  obtaining  the  discretized  flow  field,  single-  and  two-phase,  was  described  in 
Sections  1  and  2  of  this  Chapter  respectively.  To  proceed  with  numerical  conditions  based  on 
the  optimal  control  system,  it  is  assumed  that  an  entire  flow  field  Uk{xiltn)  is  known  for  all 
spatial  indices  i,  temporal  indices  n  and  variable  indices  k.  The  vast  majority  of  the  additional 
computational  complexity  required  to  satisfy  all  conditions  of  the  optimal  system  is  solving  the 
adjoint  PDE.  Since  the  method  is  nearly  identical  for  both  the  single-  and  two-phase  calculations 
it  will  be  described  first. 

A  survey  of  the  literature  on  numerical  methods  for  optimal  control  [77]-[89]  gave  insight  on 
how  to  construct  an  adjoint-based  method  of  solution.  However,  there  were  no  methods  specif¬ 
ically  for  distributed  control  with  free  initial  and  final  data  and  final  time  for  unsteady  shock 
attenuation. 

3.3.1  Discrete  Form  of  Adjoint  System 

The  continuous  form  of  the  adjoint  PDE  system  is  shown  as  the  pairing  with  SUj  in  Equation 
2.69.  The  system  is  linear  in  the  adjoint  variables.  The  time  derivative  of  the  entire  adjoint 
vector  for  all  space  at  a  discrete  time  tn  is  shown  in  column-upon-column  form  in  Equation 
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3.58.  Let  m  be  the  number  of  spatial  grid  points  and  k  be  the  number  of  adjoint  variables.  Then 
the  adjoint  vector  V  at  a  discrete  moment  in  time  will  be  of  size  km  by  1  and  the  matrices  in 
Equation  2.69  will  be  of  size  km  by  km.  A  single  component  of  the  adjoint  vector,  eg.  Vj  (x,t), 
will  be  size  m  by  1  at  each  time  step. 


( 


(  v^x)  \ 


\ 


Vifat") 


\  Vi(xm,tn)  j 
(  v2(Xl ,tn)  \ 
V2(x2,tn ) 

V  v2 (xm,tn)  j 


(3.58) 


(  Vk{Xl,tn)  \ 
Vk(x2 ,tn) 


\  V  Vk(xm,tn)  J  J 


All  of  the  matrices  in  Equation  2.69  have  a  block  diagonal  structure.  The  Jacobian  matrix  for 
all  of  space  at  tn  is  shown  in  Equation  3.59. 
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Akv(Un)  -f-)- 


/ 

^  Au(xi,tn) 

<  Alk(x1,tn) 

> 

\ 

{ 

Au{xm,tn )  j 

{ 

Alk(xm,tn)  ) 

(  Akl(xi,tn) 

(  Akk(xi,tn) 

V 

Akl(xm,tn)  j 

K 

-^4 kk  (*^rri7  t  )  / 

) 

(3.59) 

The  adjoint  PDE  is  given  in  discrete  form  with  an  explicit  integration  in  Equations  3.60  and 


3.61. 


(yrl  -  k) +R(un)  v"  =  o 


(3.60) 


R  (U)  =  A  ■  DU  + 


dA 

dx 


dA  OU  dS 
dU  dx  dU 


(3.61) 


The  matrix  D  is  made  up  of  discrete  spatial  derivative  block  matrices,  central  differencing  in 
the  domain  interior,  upwind  differencing  at  the  outlet  and  downwind  at  the  inlet.  A  single  block 
is  shown  in  Equation  3.62. 


DmxM  — 


1 

2Ax 


<  -2 
-1 


2 

0 


V 


1 

-1  0 
2 


\ 


1 

-2/ 


(3.62) 


Each  time  step  of  the  adjoint  solution  has  four  parts.  Before  time  integration,  the  matrix  R  must 
be  assembled.  Some  of  the  matrices  which  make  up  R  are  known  in  closed  form  and  require  no 
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discretized  derivative.  The  3-component  system  requires  assembling  A  and  from  the  known 
state  data.  The  7-component  system  requires,  in  addition,  assembling  which  has  terms  with 
mixed  closed  form  and  discrete  derivatives.  The  second  part  of  the  solution  requires  assembling 
the  matrix  and  vectors  that  have  discrete  derivatives  Af-  and  These  two  parts  can  be 
done  in  parallel.  The  third  part,  calculating  R,  requires  sharing  memory  between  processes 
and  does  not  lend  itself  well  to  parallelization.  With  careful  direction  of  memory,  there  is  more 
potential  in  the  assembly  of  R  for  speed  optimization  than  will  be  shown  in  the  results  in  Section 
4.4.  The  final  part  of  the  adjoint  time  step  is  the  time  integration  which  boils  down  to  matrix 
addition  and  matrix-vector  multiplication  for  the  explicit  scheme.  These  operations  are  known 
to  be  adaptable  to  parallelization  in  a  straightforward  way. 

For  adjoint  calculations  of  a  scalar  PDE  with  a  discontinuity,  it  has  been  shown  [91]  that  a 
relaxed  system  with  second  order  dissipation  will  recover  the  non-linear  PDE  in  the  limit  of 
vanishing  viscosity.  A  small  numerical  viscosity  can  stabilize  the  adjoint  solution.  These  ideas 
have  been  extended  to  fluid  dynamics  systems  [92]  and  are  implemented  in  the  current  work, 
in  a  manner  which  maintains  consistency,  for  both  the  single-  and  two-phase  numerical  adjoint 
solutions. 

3.3.2  Formulation  and  Solution  of  Single-Phase  Control  Problem 

The  goal  of  the  calculation  is  to  minimize  the  finite  approximation  to  the  cost  functional  J  while 
also  converging  to  a  solution  that  satisfies  all  of  the  necessary  conditions  of  the  optimal  system. 
The  cost  functional  from  Equation  2.79  is  given  in  discrete  form  in  Equation  3.63. 

J  =  E^1E^=1|^(a;i,fn)2AxAf  +  E Xi£^max  (Pg(xi,tN)  -  Q(xi),  o)  Ax  (3.63) 

Henceforth,  i  will  replace  j  as  the  spatial  index.  The  state  vector  U  is  constrained  to  the  Euler 
System  with  a  distributed  control  z  in  the  energy  balance  equation. 

U  =  (p(xh  tn),pu(xi ,  tn),pE(xi,  tn))T  (3.64) 
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f  P  \  (  Pu  \  (  0  > 

ift  pu  +~L  pu2+p  =  °  ^3-65) 

\pE  ]  \u(pE  +  P)  )  \z(xhtn)  j 

2 

pE(Xiltn )  =  pe(Xi,tn)  +  ^-( Xi,tn )  (3.66) 

P(Xi,tn)=pe(Xi,tn){  7-1)  (3-67) 

Initial  Conditions  are  stationary  air  at  sea-level  throughout  the  domain. 

p(Xi,t° )  =  1  kg/m? 

u(Xi,t° )  =  0  m/s 

P(Xi,t0)  =  10 5  Pa 

(3.68) 

pu(cci,t° )  =  0  kg/m2  ■  s 

pe(Xi,t°)  =  250000  J 

pE(Xi,t°)  =  250000  J 

The  control  solution  procedure  requires  a  solution  of  the  flow  under  the  influence  of  the  Ith 
control  iterate  zl.  Therefore,  the  initial  control  iterate  will  be  zero  z°  =  0.  The  first  iterate  of 
the  final  time  T°  is  chosen  so  that  the  shock  wave  is  allowed  enough  time  to  almost  traverse 
the  domain  Q.  In  this  way,  the  greatest  potential  for  the  effect  of  control  action  on  pressure 
attenuation  is  possible.  At  the  final  time  Tl  the  necessary  condition  on  the  adjoint  vector  was 
given  in  Equation  2.97  and  is  shown  in  discrete  form  in  Equation  3.69  for  all  Xi  e  Qs. 
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=  b  ■  max  (JyP  (xi,Tl^j 

-Q(xi )) 

\  U  (xi,Tl) 
•°)  2 

(7- 

-1) 

{xi,Tl) 

=  —b  ■  max  (  (p  ( Xi ,  Tl 

)  -  Q(xi) 

)  ,o)  u  (xi,Tl 

)  (7- 

-1) 

(3.69) 

(xi,Tl) 

=  b  ■  max  (Jp  (xi,Tl^j 

-Q(xt )) 

i— 1 

1 

/"o^' 

The  adjoint  PDE  is  then  solved  backward  in  time  as  described  in  Section  3.3.1.  This  gives  all 
of  the  adjoint  variables  Vj. (xj,  tn )  over  time  on  the  same  discrete  grid  as  the  flow  variables.  A 
new  control  iterate  is  obtained  from  the  discrete  form  of  Equation  2.99. 

Szl+1(xi,tn)  =  —Vl{xi,tn)  /  a 

zl+1(xi,tn)  =  zl(xi,tn )  +  5zl+1(xi,tn) 

Physical  control  constraints  are  then  imposed.  Since  the  calculation  concerns  shock  attenuation, 
extracting  energy  from  the  gas  using  a  sink  is  of  interest.  Consequently,  zl  <  0  for  all  space  and 
time  is  enforced.  It  has  been  shown  that  the  adjoint  variables  will  travel  along  the  characteristics 
of  the  flow  in  the  opposite  direction.  This  means  that  the  calculation  will  suggest  putting  control 
ahead  of  the  shock  wave  which  we  are  not  interested  in  since  droplet-shock  interactions  only 
take  energy  out  of  the  gas  behind  the  shock  wave.  Therefore,  Equation  3.72  gives  the  restrictions 
on  the  control. 


zl+1(xi,tr‘ 


0  if  zl+1(xi,  tn)  >  0 
0  if  Xi>  S  ■  tn 


(3.71) 


Lastly,  the  final  time  is  updated  according  to  Equation  3.72  and  3.73.  It  is  important  to  mention 
that  changing  the  final  time  does  not  change  the  number  of  time  steps  N.  This  means  that  the 
time  step  will  be  slightly  different  between  solution  iterates  but  will  still  obey  the  CFL  condition. 
This  also  means  that  control  strengths  defined  on  the  time  discretization  of  a  previous  iterate 
will  be  shifted  when  applied  to  the  current  flow  solution.  Overall  convergence  of  the  algorithm 
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ensures  that  these  differences  are  not  detrimental  and,  as  the  final  time  converges,  so  will  these 
discrepancies. 


f(Tl)  =  mia-z(xi,tNfAx+ 


( 


b(Pg(xi,tN )  -  Q(Xi)) 


(7  - 1)  • ' 


U, 


Ug  (Xh  pg  ( Xi ,  tN^j  ~  Pg  (x{  ,  A 

2  At 

/  , N\  PgU9  fc, 

^i,rj - - - - — — - - - —\ 


N-A  \ 


At 


PgEg  (Xi,  tN )  -  PgEg  (^Xi ,  1 

At 

Pg  (Xi,  -  Pg  (xt,  t^) 

At  ^ 


]  + 


Aa; 


(3.72) 


,  zAi 

df  , 

i<T'> 


(3.73) 


As  long  as  each  iterate  of  the  final  time  Tl+l  is  not  too  far  from  the  optimal  value  T* ,  the 
variable  will  converge  in  the  optimal  control  solution  procedure  and  the  properties  in  Equation 
3.74  will  be  true.  The  definition  of  ^  and  its  discretization  are  given  in  Appendix  C. 


lim/^-oo  Tl+1  =  T* 
lim^0O/(T'+1)  =  f  (T*)  —  0 


The  discrete  form  of  the  functional  /  from  Equation  2.98  is  given  in  Equation  3.72,  using  the 
definition  of  the  derivatives  of  pressure  from  Equation  2.97.  The  overall  algorithm  which  is 
used  to  satisfy  all  of  the  above  necessary  optimal  conditions  is  shown  in  block-diagram  form  in 
Figure  3.3 
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Blast  wave  simulation  data  at  a  fixed  point  in  space  will  be  provided  as  an  inlet  boundary 
condition.  If  the  state  of  the  flow  is  supersonic,  the  condition  will  be  explicitly  applied.  If  the 
flow  is  subsonic,  a  non-reflecting  boundary  condition  of  the  outgoing  u  —  c  characteristic  wave 
will  be  enforced  [93],  [94]. 
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Figure  3.3:  Block  diagram  of  solution  procedure  for  single-phase  control  calculation 
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3.3.3  Formulation  and  Solution  of  Two-Phase  Control  Problem 


The  solution  procedure  for  the  two-phase  control  problem  is  very  similar  to  that  of  the  single¬ 
phase  control  problem.  The  major  distinction  is  that  an  adjoint  PDE  must  be  solved  backward 
in  time  to  fix  an  initial  condition  on  a  free  state  variable  rather  than  to  calculate  a  distributed 
control  over  time.  The  cost  functional  J ,  given  in  Equation  2.100  and  shown  in  discrete  form 
in  Equation  3.76,  will  be  minimized  while  all  constraints  of  the  two-phase  model  are  obeyed. 


J  =  ^i^aKxi,t0)Ax  +  Sx.ef2s  (Pg(xi,tN)  -  Q{xi),  o)2  Ax  (3.75) 

There  is  no  running  penalty  necessary  since  the  cost  associated  with  control  action  is  penalized 
at  the  initial  time  and  the  cost  associated  with  missing  the  target  pressure  profile  need  only  be  as¬ 
sessed  at  the  final  time.  Recall  the  state  vector  in  primitive  form,  U  =  (ag,  pg,  ug,  Pg,  pi,ui,  P{)T 
and  the  constraint  on  the  free  initial  data  ag  +  ai  =  1  and  0  <  a:/  <  1,  throughout  the  simulation 
domain  Q.  The  rest  of  the  state  variables  have  fixed  initial  conditions  that  are  given  in  Equation 
2.104.  The  dynamical  constraints  of  the  two-phase  system  are  solved  in  their  conservative  form 
as  in  Equation  2.24,  but  since  it  is  convenient  to  define  the  target  state  in  terms  of  the  pres¬ 
sure,  and  the  Jacobian  has  already  been  provided  [31],  the  resulting  necessary  conditions  are 
simplified.  The  adjoint  PDE  is  derived  from  the  Jacobian  form  in  the  primitive  basis  shown  in 
Equation  2.51.  The  non-trivial  elements  of  the  matrices  are  given  in  Appendices  A  and  B. 

Analogous  to  Equation  3.69  the  necessary  condition  on  the  adjoint  vector  at  the  final  time  is 
given  for  all  x.t  G  Vts. 


b  ■  max  ^P  (xi,Tl^j  —  Q(xi)^j  ,  o)  k  —  4 
0  k  ±  4 


(3.76) 


With  the  final  time  condition,  the  adjoint  PDE  can  be  solved  backward  in  time  to  arrive  at  a 
value  for  Vj(xi,  0).  From  the  necessary  condition  given  in  Equation  2.107,  the  volume  fraction 
of  liquid  at  each  point  in  space  will  be  iterated,  within  the  constraints,  via  a  Newton  method 
shown  in  Equation  3.77. 
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(3.77) 


To  solve  for  the  optimal  final  time,  the  Transversality  Condition  for  the  two-phase  system,  given 
in  Equation  2. 1 15,  is  define  as  a  continuous  function  f(T)  with  the  final  time  as  the  independent 
variable.  Then  f(T*)  =  0  at  the  optimal  final  time  T* .  As  in  the  single-phase  case,  T*  can  be 
solved  iteratively  with  the  Newton  method  from  Equation  3.72  using  the  discretization  of  /  in 
Equation  3.78. 


,  ,  ,  \  ,  P„  (xi,  tN)  —  P„  (xi,  t1 

f{T)  =  (pg  (n,tN)  -  q(Xi))  — L^y — 


(3.78) 


The  definition  of  ^  and  its  discretization  are  also  given  in  Appendix  C.  The  overall  algorithm 
which  is  used  to  satisfy  all  of  the  above  necessary  optimal  conditions  is  shown  in  block-diagram 
form  in  Figure  3.4. 
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Figure  3.4:  Block  diagram  of  solution  procedure  for  two-phase  control  calculation 
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CHAPTER  4: 
Results  and  Discussion 

4.1  IOP  Simulations  in  2D  using  FASTRAN 


Figure  4.1:  Simulated  Shuttle  IOP:  Mach  number  at  1.2  ms 


3  3SW  <12 < 


Figure  4.2:  Simulated  Shuttle  IOP:  Mach  number  at  4  ms 


J.55K-02I 


Figure  4.3:  Simulated  Shuttle  IOP:  Mach  number  at  10  ms 

Data  on  the  Shuttle’s  grain  and  chamber  pressure  [95]  was  given  as  input  to  Cequel,  [96]  a  code 
for  steady  state  rocket  property  calculations.  The  output  gives  the  temperature  and  gas  velocity 
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at  the  nozzle  exit  plane  for  a  given  pressure  ratio.  Initially,  ambient  conditions  inside  the  domain 
are  present.  The  ignition  sequence  was  simulated  in  2D  using  the  ESI-Fastran  Solver  [97]. 


Figure  4.4:  Simulated  Shuttle  IOP:  pressure  at  1.2  ms 


Figure  4.5:  Simulated  Shuttle  IOP:  pressure  at  4  ms 


Figure  4.6:  Simulated  Shuttle  IOP:  pressure  at  10  ms 


Constant  mass-flow  boundary  conditions  equivalent  to  the  steady  state  exit  plane  of  the  rocket 
nozzle  on  the  shuttle  were  used  in  the  bottom  center  of  the  domain  on  the  right  face  of  the  step 
as  shown  in  Figure  4.1.  Mach  number  is  depicted  in  three  snapshots  in  Figures  4.1,  4.2  and  4.3 
and  pressure  in  Figures  4.4,  4.5,  and  4.6.  The  last  frame  is  roughly  10  ms  after  ignition.  The 
bottom  left  edge  of  the  domain  represents  the  rocket  body  while  the  bottom  right  edge  is  the 


70 


centerline  of  the  normal  to  the  nozzle  exit  plane  and  a  symmetry  boundary.  All  other  edges  are 
non-reflecting  boundaries. 

Flow  conditions  over  time  were  recorded  at  two  locations  marked  in  Figure  4.1.  Point  1  is 
near  the  rocket  body  2.5  meters  above  the  nozzle  and  Point  2  is  1.5  meters  along  the  symmetry 
boundary  and  the  plane  normal  to  the  nozzle  exit.  The  conditions  at  the  recorded  locations  are 
used  as  the  boundary  conditions  in  both  the  single-  and  two-phase  control  calculations. 

Figure  4.7  shows  the  flow  conditions  over  time  at  Monitor  Point  1  near  the  rocket  body  and 
Figure  4.8  shows  the  flow  conditions  over  time  for  Monitor  Point  2  directly  downstream  of  the 
nozzle. 
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Figure  4.7:  Flow  conditions  over  time  for  Monitor  Point  1 


It  is  worth  mentioning  that  the  two-dimensional  flow  cannot  be  truly  replicated  by  an  inlet 
boundary  condition  to  a  single-dimensional  calculation.  Neglecting  motion  in  the  transverse  di- 
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Figure  4.8:  Flow  conditions  over  time  for  Monitor  Point  2 

rection  while  taking  (p,  u,  P )  explicitly  means  that  the  resulting  driver  gas  has  a  lower  temper¬ 
ature.  This  is  not  discouraging  since  intuition  would  suggest  that  hotter  gas  has  more  potential 
to  vaporize  water  droplets  and  hence  by  the  nature  of  the  control  action,  there  is  more  potential 
for  effectiveness.  The  goal  is  to  be  able  to  handle  a  range  of  blast  waves  that  will  plausibly  be 
seen  in  an  IOP  environment. 
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4.2  Single-Phase  Control  Calculation  Results 


Space  (m)  Time(s) 


Figure  4.9:  Monitor  Point  2,  Pressure  profile  at  final  time,  no  control  used. 

These  results  were  published  by  the  author  [98].  The  solution  procedure  in  Figure  3.3  assumes 
a  given  target  state  Q(x)  at  the  final  time  T*.  To  illustrate  trends  in  the  optimal  control  solution, 
prescribing  a  consistent  and  meaningful  target  state  or  sequence  of  target  states  is  a  necessity. 
Figure  4.9  shows  the  uncontrolled  pressure  profile  over  time.  An  example  of  a  target  state, 
shown  in  red  in  Figure  4.10,  has  85%  of  the  amplitude  of  the  final  pressure  profile  when  no 
control  is  used,  OPq,  shown  in  blue.  This  attenuated  pressure  profile  replaces  the  uncontrolled 
pressure  profile  from  just  behind  the  shock  front  to  30  cm  from  the  inlet  boundary.  The  tar¬ 
get  is  then  linearly  increased  to  the  magnitude  of  the  uncontrolled  pressure  profile  over  10cm. 
Upstream  of  the  shock,  the  target  state  is  such  that  only  pressures  greater  than  the  target  pres¬ 
sure  are  penalized.  This  will  define  Qs  in  the  final  time  penalty.  The  nonlinear  nature  of  shock 
waves  means  that  not  all  target  states  may  be  possible  for  the  given  initial  conditions  and  bound¬ 
ary  conditions.  In  addition,  the  purpose  of  the  calculation  is  to  calculate  control  solutions  which 
decrease  overpressure  and  therefore  the  cost  functional  should  not  penalize  pressures  below  the 
target  pressure.  Furthermore,  error  near  the  shock  front  is  not  directly  penalized  but  is  still 
minimized  via  the  iterative  update  to  the  final  time.  This  avoids  taking  a  variation  of  the  state 
variables  near  discontinuities,  which  will  violate  the  assumption  of  a  small  variation.  In  each 
of  the  given  examples  the  weighting  constants  a  and  b  are  10-6  and  104  respectively.  Since  b  is 
much  larger  than  a,  minimizing  J  is  dominated  by  minimizing  the  final  time  penalty.  The  value 
for  b  is  such  that  the  magnitude  of  the  third  adjoint  variable  V%(x,  t ),  and  thus  the  control  action 
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z(x,  t ),  are  sufficient  to  cause  noticeable  attenuation  to  the  blast  wave. 
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Figure  4.10:  Monitor  Point  2,  Optimal  pressure  profile  at  final  time,  target  state  Q(x)  ~ 
85%  OP0 


Time  (s) 


2  0 

Space (m) 


2 


Figure  4.11:  Monitor  Point  2,  Optimal  distributed  control, Q(x)  ~  85%  OP0 

The  desired  state  cannot  be  directly  set  as  described  but  must  be  gradually  scaled  down  to  the 
desired  overpressure  to  assure  convergence  of  the  algorithm.  It  was  found  that  when  starting 
with  an  overpressure  of  8:1  from  Monitor  Point  2  data,  setting  a  target  state  with  85%  of  the 
uncontrolled  overpressure,  as  in  Figures  4.10,  was  the  most  attenuation  to  the  shock  that  could 
be  desired  and  still  maintain  convergence  of  the  the  algorithm.  The  reason  for  this  can  be 
understood  in  the  nature  of  the  optimal  system.  The  optimal  conditions  are  all  coupled  and 
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Figure  4.12:  Monitor  Point  2,  Optimal  pressure  at  final  time,  Q(x)  ~  70%  OP0 
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Figure  4.13:  Monitor  Point  2,  Optimal  distributed  control,  Q(x)  ~  70%  OP0 


inter-dependent  such  that  a  calculated  blast  wave  under  no  control  action  is  in  no  way  close  to 
the  optimal  pressure  profile  of  a  blast  wave  under  considerable  damping  control  action  so  as  to 
render  a  blast  wave  that  matches  the  target  state  with  a  significantly  diminished  overpressure. 
When  the  final  time  solution  is  not  nearly  optimal,  then  first-order  corrections  to  these  terms 
are  not  enough  to  exploit  in  an  iterative  sense  toward  finding  an  optimal  control  satisfying  all 
necessary  conditions. 

Figure  4.11.  is  the  corresponding  distributed  optimal  control  solution.  Figures  4.12  and  4.13 
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are  the  optimal  pressure  and  control  solutions  respectively  for  a  target  state  that  has  70%  of  the 
uncontrolled  overpressure.  Figures  4.14  and  4.15  are  the  optimal  pressure  and  control  solutions 
respectively  for  a  target  state  that  has  50%  of  the  uncontrolled  overpressure. 
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Figure  4.14:  Monitor  Point  2,  Optimal  pressure  at  final  time,  Q(x)  ~  50%  OPq 
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Figure  4.15:  Monitor  Point  2,  Optimal  distributed  control,  Q(x)  ~  50%  OP0 

The  top  plot  in  Figure  4.16  is  the  logarithm  of  the  cost  functional  J  over  solution  iterations. 
The  jumps  in  J  every  25th  iteration  indicate  that  the  target  state  Q(x)  has  been  redefined  with 
a  lower  overpressure.  The  bottom  plot  in  Figure  4.16  shows  the  corresponding  final  time  T* 
iterates  with  T°  =  2  ms. 
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A  negative  control  in  the  energy  equation  can  be  thought  of  as  an  internal  energy  sink.  Tem¬ 
perature  is  directly  related  to  internal  energy,  so  the  control  action  cools  the  gas.  A  cooler  gas 
has  a  lower  speed  of  sound.  Consequently,  it  will  take  longer  for  the  shock  front  in  the  wave 
interacting  with  control  to  reach  the  same  point  in  space  that  the  shock  without  control  reached. 
Additionally,  from  a  hyperbolic  wave  theory  point-of-view,  a  non-linear  wave  travels  slower  if 
it  has  less  amplitude.  The  control  action’s  purpose  is  to  attenuate  amplitude  and  this  necessar¬ 
ily  slows  the  wave  down.  The  Transversality  Condition  is  very  sensitive  to  error  at  the  shock 
front  since  this  is  where  the  time  rate  of  change  of  the  pressure  is  greatest.  It  would  be  a  diffi¬ 
cult  matter  to  get  significant  attenuation  in  a  converged  solution  with  a  fixed  final  time  for  this 
reason. 


Figure  4. 16:  Monitor  Point  2,  Cost  functional  (top)  and  final  time  (bottom)  vs.  iteration,  Q(x)  ~ 
50%  OP0 

In  each  example,  a  physical  restriction  on  the  control  is  imposed  such  that  energy  can  only  be 
taken  out  of  the  gas  behind  the  shock  wave.  This  more  accurately  portrays  how  discrete  droplet 
sprays  will  sink  energy  from  the  gas  via  vaporization.  The  plots  shown  in  Figures  4.17  and  4.18 
are  distributions  of  the  energy  equivalent  vaporized  water  mass  to  the  optimal  control  solutions. 
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By  dimensional  analysis,  inspection  of  the  energy  balance  equation  indicates  that  the  units  of 
z*(x,  t )  are  Watts.  Integrating  the  control  solution  from  [0,  T*]  gives  an  energy  distribution  in 
space.  This  energy  can  be  directly  equated  to  an  amount  of  water  mass  that  must  be  vaporized 
using  the  latent  heat  of  vaporization  of  water  Lhv  =  —2.26  ■  106  J /kg  at  100°  C.  Equation  4.1 
relates  the  optimal  control  solution,  distributed  in  space  and  time,  to  an  equivalent  distribution 
of  water  mass  vaporized  over  the  same  time  interval. 

1  rT* 

mH2  qv(x)  =  - —  /  Z*(x,t)dt  (4.1) 

L>hv  JU 


Figure  4.17:  Monitor  Point  2,  mH20vape(x),  Energy  equivalent  distribution  of  water  mass  vapor¬ 
ized 
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Figure  4.18:  Monitor  Point  1,  rriH2ovape{x),  Energy  equivalent  distribution  of  water  mass  vapor¬ 
ized 
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Figures  4. 17  and  4. 18  show  how  the  equivalent  vaporized  water  mass  distribution  increases  with 
increasing  shock  front  attenuation.  It  can  clearly  be  seen  that  more  water  is  required  to  attenuate 
the  stronger  shock  front  using  Monitor  Point  2  data  from  directly  in  front  of  the  rocket  nozzle 
than  from  the  weaker  shock,  Monitor  Point  1  data,  taken  near  the  rocket  body.  Water  takes  time 
to  vaporize  and  therefore  cannot  sink  energy  from  the  gas  instantaneously  and  arbitrarily  close 
to  the  shock  front.  In  practice,  water  will  take  energy  out  of  the  gas  at  least  a  meter  behind 
the  shock  front.  That  less  energetic  gas  will  then  drive  the  shock  with  less  energy  and  the 
pressure  profile  will  flatten  out  to  a  diminished  overpressure.  With  more  detailed  physically- 
based  restrictions  on  the  control  action,  this  single-phase  shock  wave  control  formulation  can 
be  made  to  better  replicate  a  two-phase  shock-droplet  interaction  control  calculation  without 
the  considerable  added  complexity  of  a  two-phase  compressible  flow  model.  Regardless,  trends 
in  the  data  seem  heuristically  correct  and  therefore  empirical  scaling  laws  in  magnitude  and 
location  may  be  sufficient. 

4.3  Two-Phase  Calculation  Results 

The  results  for  the  two-phase  control  problem  were  obtained  with  the  algorithm  shown  in  Figure 
3.4.  Figure  4.19  shows  the  optimal  initial  condition  for  the  free  water  volume  fraction  variable. 
The  curves  give  the  optimal  water  volume  fraction  distribution  for  a  given  target  state.  The 
legend  tells  what  percentage  of  the  absolute  overpressure  of  the  uncontrolled  blast  wave  was 
used  to  define  the  target  state.  In  this  case  there  is  no  artificial  upper  constraint  on  the  maximum 
value  of  the  volume  fraction,  aside  from  the  physical  limit  of  1 . 

The  limit  of  controllability  over  two  meters  is  close  to  the  target  state  76%  OP$.  For  each 
optimal  initial  water  distribution,  the  resulting  pressure  profiles  at  the  optimal  final  time  T*  are 
shown  in  Figure  4.20. 
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Figure  4.19:  Optimal  water  volume  fraction  distributions  a*(x,  0),  unconstrained  at  initial  mo¬ 
ment  in  time 
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Figure  4.20:  Optimal  final  time  pressure  profiles  resulting  from  initial  data  in  Figure  4.19 
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Figure  4.21:  Cost  functional  J  in  red  and  initial  penalty  /  over  iterations  of  solution  procedure 
in  Figure  3.4 
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Figure  4.22:  Final  time,  in  seconds,  vs  iteration  number  for  the  unconstrained  case.  The  regions 
of  greatest  slope  are  after  iterations  where  the  target  state  has  been  changed.  The  slope  tending 
toward  zero  means  that  the  solution  is  converging. 
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The  convergence  of  the  algorithm  is  shown  in  Figure  4.21.  Every  25  iterations,  the  target  state’s 
amplitude  is  ramped  down  so  that  optimal  solutions  that  yield  ever  increasing  attenuation  levels 
can  be  determined  in  a  single  execution  of  code.  The  history  of  the  final  time  variable  over  the 
iterative  solution  procedure  is  shown  in  Figure  4.22. 


Figure  4.23:  Optimal  initial  water  volume  fraction  distributions,  cei(x,  0)  <  .7 


Another  case  was  run  with  an  upper  bound  on  the  volume  fraction  of  water  at  70%.  The  calcu¬ 
lated  optimal  water  distributions  are  shown  in  Figure  4.22  with  the  resulting  pressure  profiles 
shown  in  Figure  4.23. 
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Figure  4.24:  The  optimal  pressure  profiles  at  the  final  time  resulting  from  the  control  initial  data 
in  Figure  4.23 
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Figure  4.25:  Optimal  initial  water  volume  fraction  distributions,  ai(x,  0)  <  .5 

A  final  case  was  run,  this  time  with  an  upper  bound  on  the  volume  fraction  of  water  at  50%  .The 
calculated  optimal  water  distributions  are  shown  in  Figure  4.25  with  the  resulting  pressure  pro¬ 
files  shown  in  Figure  4.26. 
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Figure  4.26:  Optimal  final  time  pressure  profiles  resulting  from  initial  data  in  Figure  4.25 
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A  few  more  isolated  results  were  obtained  in  the  interest  of  categorizing  the  significance  of 
model  parameters.  Figure  4.27  shows  the  effect  of  droplet  size  on  overpressure  attenuation. 
Each  pressure  profile  is  plotted  after  two  milliseconds  of  shock  propagation  from  the  left  bound¬ 
ary  toward  the  right.  For  a  constant  amount  of  water  mass,  more  surface  area  of  the  droplets 
are  exposed  to  the  flow  when  the  droplets  are  smaller.  This  is  verified  in  Equation  4.2  for 
mono-dispersed  droplets  with  diameter  Dp. 


Sn  —  47T  (  ^r-  )  '  N, 


Np  cx[- 


2 

Aa:3 


(4.2) 


This  result  agrees  with  empirical  results  that  more  exposed  droplet  surface  area,  the  greater  the 
overpressure  attenuation  [12].  Since  surface  area  is  maximized  with  infinitely  small  droplets, 
initial  droplet  size  isn’t  an  optimizeable  parameter.  Droplets  with  a  diameter  Dp  =  25/rm  were 
used  in  all  of  the  preceding  results  since  that  size  is  small  enough  to  have  a  significant  effect 
over  two  meters,  is  a  reasonable  size  based  on  injector  atomizer  specifications  and  prevents  the 
droplets  from  being  completely  vaporized  within  the  simulation  time  T*.  Avoiding  complete 
vaporization  is  desirable  since  a  non-zero  volume  fraction  of  either  phase  is  a  prerequisite  to  the 
model  being  used  [31]. 
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Figure  4.27:  Effect  of  variable  droplet  size 
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Water  vapor  density  was  shown  to  have  a  drastic  effect  on  results  for  optimal  water  distributions. 
According  to  the  steam  table  data  [42],  water  vapor  density  increases  by  an  order  of  magnitude 
from  ambient  pressure  to  the  8  atm  level  behind  the  shock  front.  Therefore,  if  water  vaporizes 
under  high  pressures,  much  more  mass  changes  phase  and  consequently,  the  dissipative  effect 
of  vaporization  is  much  more  significant.  The  results  in  Figure  4.28  show  that  more  than  twice 
as  much  water  is  needed  to  get  5%  attenuation  in  pressure  depending  on  what  the  gas  pressure 
surrounding  the  droplets  when  they  vaporize. 


Figure  4.28:  Effect  of  variable  water  vapor  density  due  to  gas  pressure 
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Figure  4.29:  Effect  of  maximum  of  Gaussian  water  volume  fraction  distribution  on  IOP  strength 
after  1  ms. 

Figure  4.29  shows  that  for  an  initial  water  volume  fraction  distribution  with  a  Gaussian  shape,  as 
the  water  volume  fraction  at  the  inlet  approaches  1,  the  effect  of  plume  confinement  [29]  takes 
over  and  can  actually  increase  the  jump  in  pressure  at  the  shock  front.  The  Gaussian  distribution 
is  such  that  the  maximum  exists  at  the  inlet  boundary  and  about  1%  of  the  maximum  half-way 
into  the  domain. 
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4.4  Serial  and  Parallel  Code  Run-Time  Optimization 


The  last  section  of  the  results  shows  the  run-time  optimization  performed  within  Matlab.  For 
both  the  single -phase  (3  adjoint)  and  the  two-phase  (7  adjoint)  systems,  5  cases  of  increasing 
size  were  carried  out  with  serial  and  parallel  syntaxes.  Jacket  [99]  was  used  to  run  the  3  adjoint 
calculation  in  parallel  on  Nvidia  GPUs.  Results  for  the  time  duration  needed  to  solve  the  system 
in  [T*,  0]  are  shown  in  Figure  4.30  and  those  of  the  7  adjoint  system  are  shown  in  Figure  4.31. 


Case  1: 
Case  2: 
Case  3: 
Case  4: 
Case  5: 


m  = 
m  = 
m  = 
m  = 
m  = 


100,  timesteps  =  1000 
200,  timesteps  =  2000 
300,  timesteps  =  3000 
400,  timesteps  =  4000 
500,  timesteps  =  5000 


105  discrete  points  in  space-time 
4  •  105  discrete  points  in  space-time 
9  •  105  discrete  points  in  space-time 
1.6  •  106  discrete  points  in  space-time 
2.5  •  106  discrete  points  in  space-time 


Figure  4.30:  Time  duration  for  solving  the  3-component  adjoint  system  of  the  Euler  equations 
over  the  time  interval  [T*,  0] 
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Figure  4.31:  Time  duration  for  solving  7-component  adjoint  system  of  two-phase  model  over 
the  time  interval  [T* ,  0] 
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CHAPTER  5: 
Conclusions 


A  new  iterative  solution  procedure  was  developed,  which  can  calculate  optimal  distributed  con¬ 
trol  solutions  for  systems  of  quasi-linear  hyperbolic  partial  differential  equations  with  free  initial 
and  final  data  and  final  time.  This  procedure  has  been  successfully  applied  to  single-  and  two- 
phase  compressible  gas  dynamics  in  one  dimension  with  the  goal  of  diminishing  overpressure 
at  the  shock  front  of  a  blast  wave  generated  by  an  ignition  overpressure.  Examples  of  opti¬ 
mal  attenuation  to  blast  waves  typically  encountered  in  the  launch  environment  of  the  Shuttle’s 
SRBs  during  an  ignition  are  given.  The  single-phase  control  solutions  can  be  seen  as  a  spatial 
distribution  of  energy  equivalent  vaporized  water  mass  required  to  achieve  a  given  level  of  over¬ 
pressure  reduction.  Results  for  these  mass  distributions  are  shown  for  inlet  boundary  conditions 
like  those  given  in  the  monitor  point  data. 

The  same  methodology  for  solving  the  control  system  is  implemented  on  a  two-phase  model. 
The  control  action  takes  the  form  of  the  free  initial  data  describing  a  distribution  of  water 
droplets.  Optimal  water  volume  fractions  are  calculated  for  increasing  levels  of  attenuation. 
Cases  where  the  maximum  volume  fraction  of  water  is  restricted  to  50%  and  70%  are  also 
presented. 

The  results  give  several  key  insights  relevant  to  implementing  water  injection  systems  for  blast 
wave  attenuation.  Smaller  droplets  will  vaporize  quicker  since  the  total  surface  area  exposed 
to  the  flow  is  larger  [12],  resulting  in  cooler  gas  driving  a  shock  with  an  attenuated  jump  at 
the  shock  front.  This  would  suggest  that  large  regions  of  space  completely  filled  with  water 
are  sub-optimal  for  blast  attenuation  since  connected  streams  or  ligaments  of  water  expose  less 
surface  area  to  the  gas. 

The  degree  of  attenuation  depends  largely  on  the  rate  of  mass  changing  phase  from  liquid  water 
to  gaseous  water  vapor  due  to  forced  vaporization  m,  which  is  equal  to  the  product  of  the  water 
vapor  density  and  the  rate  of  change  of  the  volume  occupied  by  the  water  droplets.  In  the  high 
pressure  and  high  temperature  region  behind  the  shock  front,  water  vapor  density  is  high  and  the 
rate  of  vaporization  is  high  as  well  which  both  contribute  to  a  large  value  for  fa.  As  the  control 
takes  effect,  the  pressure  and  temperature  of  the  gas  both  decrease  meaning  that  the  effects  are 
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being  felt  by  the  shock  front  at  a  slower  rate.  These  results  show  that,  in  general,  stronger 
shocks  can  be  attenuated  more  rapidly  than  weaker  shocks  via  water  droplet  vaporization. 

Similar  to  work  that  has  identified  plume  confinement  as  being  disadvantageous  to  blast  attenu¬ 
ation  [29],  the  results  also  show  that  when  a  shock  propagates  through  a  water  volume  fraction 
very  close  to  1  over  a  sufficient  interval  in  space,  so  much  liquid  water  is  converted  to  gaseous 
water  vapor  that  the  increase  in  gas  density  becomes  more  significant  to  the  pressure  than  the 
cooling  of  the  gas  due  to  vaporization.  This  shows  that  for  a  given  interval  of  space  there  is 
a  limit  to  the  controllability.  In  a  two  meter  domain,  76%  OP0  was  the  target  with  the  low¬ 
est  overpressure  for  which  the  solution  still  converged.  The  case  with  the  volume  fraction  of 
water  restricted  to  a  maximum  of  ai  <  70%  converged  with  the  lowest  target  overpressure 
of  80%  Or\).  With  the  water  restricted  to  below  oti  <  50%  by  volume,  the  most  attenuation 
achievable  with  a  two  meter  domain  is  82.5%  OP0. 

The  code  developed  by  the  author  is  still  in  the  prototype  stage,  however,  the  run-time  opti¬ 
mization  results  show  that  certain  parts  of  the  code  have  been  optimized  within  Matlab  script¬ 
ing.  Implementation  of  the  two-phase  calculation  on  parallel  GPUs  using  Jacket  or  through 
other  GPU  platforms  would  make  the  code  a  feasible  tool  to  use  in  IOP  and  other  two-phase 
unsteady  shock  wave  attenuation  calculations.  The  framework  of  the  control  formulation  and 
method  of  solution  will  remain  unchanged  while  allowing  for  added  degrees  of  freedom  like 
poly-dispersed  droplets,  droplet  breakup,  surface  tension  and  multiple  dimensions. 
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APPENDIX  A: 

Calculation  of  the  elements  of  for  the 

two-phase  adjoint  system 


Row  1: 


Row  2: 


Row  3: 


d  , 

o — 
OUk 

=  m 

dx 


dV^do^  dV^dpg  dV^du^  |  dVjdpi  |  dVjdui 

dag  dx  dpg  dx  dug  dx  dpi  dx  dui  dx 


(A.l) 


d  _  duk 
duk  21  dx 


Pg^  (dVi,  _  ChiA  Vj-Ug  /  dpg^ 
ag  \dx  dx  )  dg  y  dx 


Pa  daa 
dg  dx 


(A. 2) 


d  duk  _  dug 

77 - ^22  '  “Tj —  —  “Tj — 

OUk  ox  ox 


(A. 3) 


d  duk  _  dpg 

o - ^23  * 

OUk  ox  ox 


(A.4) 


d  _  duk 

duk  31  dx 


-1  ^  +  (P  _p.)  f  1  fr*g  ,  1  dpg 

dgPg  dx  dx  9  1  \dg  dx  pg  dx 


(A. 5) 


d  duk  dug 

~r\ - A33  '  — -  —  — — 

OUk  ox  ox 


(A. 6) 


d  _  duk  =  -Idpg 

duk  34  dx  p2  dx 


(A.7) 
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Row  4: 


*9  |  _  PgCgi  f  dVi  d)Ug  Vi  Ug  d(Xg  \  Vi  Ug  /  ,  d Pi  d)Pg 

duk  dx  ag  y  dx  dx  ag  dx  )  ag  \  '  dx  dx 


d  duk  dPg 

~r\  4l43  •  —  =  7—^ 

UUh  ox  ox 

d  duk  _  dug_ 

duk  44  dx  dx 

Row  5: 

—A  dUk  -  Pl  ( —  ®Ul  ^  |  Vi~Ul  (dpi  !  Pi  ®a9 
duk  51  dx  1  —  ag  y  dx  dx  J  1  —  ag  y  dx  1  —  ag  dx 

d  duk  _  dui 

"7 - ^55  ‘  “T-  ~  “T- 

OUk  ox  ox 

d  .  duk  _  dpi 

o - ^56  ‘  “T-  —  “7 

otik  dx  ox 


Row  6: 


d  _  duk 
duk  61  dx 


-1 

(!  -  ag)  Pl 


<m 

dx 


m 

dx 


+  (Pi  ~  Pd 


-1  dag  +  1  <9p;\ 
1  —  ag  dx  pi  dx  J 


d  duk  _  dui 

7 - ^66  ' 

OU  k  ox  ox 

d  _  duk  _  -1  dpi 
duk  67  dx  pf  dx 


(A. 8) 

(A. 9) 

(A. 10) 

(A.  11) 

(A. 12) 

(A. 13) 

(A. 14) 

(A. 15) 

(A. 16) 
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Row  7: 


_d_A  duk  _  Plcl  (DVj  dm  Vj-mdaA  v—uO  dp%  dPA 

duk  71  dx  1  —  ag  y  Dx  dx  1  —  ag  dx  J  1  —  ag  y  ^  Dx  dx  J 

(A. 17) 


d  duk  _  dPi 

^ - A76  •  — —  —  7 l~pr— 

UUk  ox  ox 


(A.  18) 


d  duk  _  d]M 
duk  7  dx  dx 


(A. 19) 
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APPENDIX  B: 

Calculation  of  the  elements  of  #  for  the 
two-phase  adjoint  system 

dS 

Solution  of  the  adjoint  system  of  PDEs  requires  deriving  the  matrix  where  S  (U)  is  the 

source  vector  for  interactions  between  the  two  phases. 


li  ( Pg  -  Pi)  +  Aag 


•mVi  +  Fd 

S ( U)  =  rn  ( Lhv  +  Ej )  +  Qi  +  FdVi  —  fiPi  ( Pg  —  Pi) 


(B.l) 


-mVi  -  Fd 

v  —m  (Lhv  +  Ei)  -  Qi  -  FdVi  +  /xP*  ( Pg  -  Pi)  j 

Before  calculating  derivatives  with  respect  to  state  variables,  it  must  be  made  clear  which  terms 
are  explicit  functions  of  which  variables.  A agvap  ( U ) ,  m  ( U ) ,  Fd  (U) ,  Qi  (U) ,  Vt  (U) ,  Pi  (U) , 
and  Ei{U)  yet  only  E,  is  an  explicit  function  of  all  seven  state  variables.  Therefore,  it  is  worth 
determining  only  the  derivatives  that  exist  ahead  of  time  to  more  precisely  define  the  elements 

J  S 

ofwr 

The  easiest  term  to  differentiate  is  the  interface  pressure. 


Pi  —  agPg  +  (1  —  Oig)  Pi 


(B.2) 


Therefore,  P,  =  P,  (ag,  P,r  P{)  and  the  derivatives  are: 


dP 

~d^n=P9~Pl 


(B.3) 
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(B.4) 


d  Pi 
dPg 


=  a 


9 


dPi 

d  Pi 


1  OLg 


(B.5) 


The  interface  velocity  is  a  function  of  a  subset  of  state  variables  as  well:  V*  =  Vt  (ag,  pg,ug,  pi,ui) 
The  derivatives  are: 


dag 


PgUg  ~  PlUl 

OL gPg  T  (1  Otg)  Pi 


OLgPgUg  T  (1  otg)  piUi  ,  % 

(cLgPg  +  (1  -  ag)  Plf  '  {P9  ^ 


(B.6) 


dVj  _  agUg _ ^gPgUg  +  (l  ~  Qg)  PlU\ 

dp g  agpg  +  (1  -  ag)  pt  (agpg  +  (1  -  ag)  ptf 


(B.7) 


dVj  _  (1  -  ag)  Ul  OLgPgUg  +  (l  ~  Qg)  PjUj  ,l_a 

dpi  agPg  +  (1  -  ag)  Pi  (agPg  +  (1  -  ag)  pif 


(B.8) 


dVi  _  OlgPg 

9Ug  CXgPg  T  (l  Ol g  )  Pi 


(B.9) 


dVi  _  (1  OLg)  Pg 

dui  Ol g P g  +  (1  -  Otg)  pt 


(B.10) 


The  drag  force  Fd  ( U )  =  Fd  (ag,  pg,  ug,  m ) 


dF d_ 
dag 


3 

4 


ddJ dPgD g  (llg 


Ul f 


(B  .11) 


dFd 

dpg 


OLg)  ( Ug  -  Uif 


(B.12) 


dFd 

dug 


\c0„Dl  (1 


OLg)  ( ug  ui) 


(B.13) 
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3 

^CdPgDp  (1  —  CXg)  ( Ug  —  III ) 


(B.14) 


dFd 

dui 


The  derivatives  of  m  and  A agvap  are  closely  related.  Recall, 


m  ( U )  =  Aagvap  ( U )  •  PH2ov(Pg)  (B.15) 

where  Ph2ov  is  the  density  of  water  vapor  and  is  only  a  function  of  gas  pressure.  It’s  derivative 
is  easily  calculated  from  the  quadratic  fit  to  the  empirical  steam  table  data  given  in  Equation 
2.38. 


dpH2°v  «  f .5368  -  .0048688  •  kg/m3  (B.16) 

dP g/ atm  \  atm) 

It  suffices  to  calculate  derivatives  of  Aagvap  with  respect  to  the  state  variables.  The  vaporization 
model  is  based  on  the  Empirical-Beta  Vaporization  Law  which  states: 


(Dp1)2  =  ( D ;)2  -  A tp  (t;) 


(B.17) 


and 


-  ((^")2  -  (ft>’  -P-'))372)  (B-1?) 

Tg  =  where  R  is  the  specific  gas  constant  for  air,  so  Aagvap  ( U )  =  Aagvap  (pg,  Pg)  .  The 
derivatives  are: 


d  A  a 


gvap 


9pg 


Npit  /3\ 
6 Ax3  V2  ) 


((D”)2-A(/?) 


1/2 


(Af) 


dp 

dpg 


(B  .19) 


and  similarly, 
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dA  Olgyap  _  Np  7T  f3\  (  f  ^n\  2 


dPn  6  Ax3  V2 


(d;)  -At f3)  (At) 


<9P„ 


(B.20) 


The  vaporization  term  is  given  by  Equation  B.21: 


P  ( Pg ,  P9)  =  p0  ^1  +  Bi  -  Trnir}j  j  /im2/,s 


(B.21) 


where  the  constants  are  given  as  P0  =  7600,  Pi  =  7.4-10  7 ,  B2  =  2.7548  and  Tmin  =  300 
degrees  Kelvin.  Therefore,  the  derivatives  are: 


r>  r>  r>  (m  m  \i?2  —  1  ( 


dpt 


—  —B0BiB2  (Tg  —  Tn 


Pa, 


(B.22) 


i~>  r>  r>  trn  m  \i?2  — 1  (^9 


d  Pn 


—  BqBiB2  {Tg  —  T„ 


Pn 


(B.23) 


The  total  energy  at  the  inerface  is  simply  defined  in  the  conservative  basis: 


Bi  —  OigE g  +  (1  —  Otg)  E[  (B.24) 

But  not  so  easily  in  the  primitive  basis.  Using  the  equations  of  state: 

=  “*  -  p)  +  d  - «»)  (^Ar)  -  p)  <B-25) 

In  this  basis,  the  interface  velocity  is  a  function  of  all  primitive  variables  and  therefore  there  are 
seven  non-zero  derivatives  to  calculate: 


dE1 

d(Xg 


Eg  -  El 


(B.26) 


114 


(B.27) 


d  Ej  _  - OLgPg 

dpg  Pg  ilg  -  !) 


d Ej  _  -  (1  -  gg)  (Pi  -  jitti) 

dpi  pj  (' 7/  -  1) 


(B.28) 


d_K 

dug 


—  agug 


(B.29) 


d_E± 

dui 


(1  -  ag)  ui 


(B.30) 


d  Ej  _  ag 

dPg  Pg  (7 g  ~  1) 


(B.31) 


dEj  _  1  -  OLg 

dPi  pi  (7/  —  1) 


(B.32) 


Lastly,  the  derivatives  of  the  rate  of  convective  heat  transfer  between  the  two  phases  at  the 
interphase  Qi  (U)  are  calculated. 


Qi  (U)  =  h  -  Sp  -  (Tg  -  Tj) 

=  ^Sp-(Tg-Ti) 

=  (2  +  .459PrV3jRe-)^.47r(t)2-(^ 


1  /  Pi+im 
CV1 


=  2  + 


45g  ^  DpPg(ug~ul ) 


1  ( Pi+im 
Cvi  vMtz-i) 


Therefore  Qi  =  Q,  (pg,  ug,Pg,pi,Pi)  and  the  derivatives  are: 


(B.33) 


dQj 

dpg 


XnDp 


Ti)  -  Nu1^] 
Pg) 


(B.34) 
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(B.35) 


d  Qj 

DUg 


XnDp 


dNu 

dug 


(T9 


Ti ) 


dQi 

dPg 


\ixDvNu 


Tg_ 

P9 


(B.36) 


dQi 

dpi 


T 

XnDpNu  — 
Pi 


(B.37) 


dQi 

dui 


XnDp 


dNu 

dui 


(T9 


Ti ) 


(B.38) 


dQi 

m 


XnDpNu— 

Pi 


T{ _ 

-  im 


And  it  can  easily  be  seen  that  the  derivatives  of  the  Nusslets’  Number  are: 


(B.39) 


^  =  .458  •  .55  •  prP3Re55-lDp^9  Ul) 

dpg  P 


(B  .40) 


^  =  .458  •  .55  •  Pr^Re  55-1^  (B.41) 

dug  p 


^  =  -.458  •  .55  •  Pr^Re  55-1^-  (B.42) 

OUi  p 

Now,  with  all  of  the  derivatives  of  each  each  individual  function  of  the  state  vector  that  compose 
the  source  vector,  the  on-zero  elements  of  the  matrix  can  be  calculated. 

Row  1: 


dSi 

dpg 


Npn  /— 3\ 
6Aa:3  V  2  ) 


{(Dr)2  -  Ai/3”) 


1/2 


d/3n 

dpg 


(B.43) 
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(B  .44) 


Row  2: 


Row  3: 


dSs  =  )hcm  m 
dag  dag  dag 

dS3  =  dm  .  dVi_  dFd 

Dpg  dpg  1  m0pg  Opg 

dS3  _  .  dV, 
dpi  m  dPl 

dSs  =  .JV,  dF^ 

dUg  dUg  dUg 

dSs  =  ihdv^  dFd 

dui  dui  dui 


(B.45) 

(B  .46) 

(B.47) 

(B.48) 

(B  .49) 

(B.50) 

(B.51) 

(B.52) 


dS3  dm 

dPg  =  dP, 


(B.53) 
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Row  4: 


dS4  dEi  <9F,  ^  dVi 

=  m-^  +  ~^Vt  +  Fd -E  -  /i  P3  -  P*  2 
ao;3  aa9  aa9  aa9 


dS4  _  dm  ^  ,  „.dEt  <9Q*  (9PdT/  <91/; 

"n —  —  "n —  \Lhv  +  Ei)  +  m— h  "7: h  t: — Vi  +  Fd— — 

Pg  ^Pg  VPg  dpg  &Pg  Vpg 


dS4  _  .  d Et  DQ i  dFd  dVi 

a  —  m  a  t~  a  h  a  V,  +  Pi  a 

CTUg  (i/Ug  OUg  OUg  OUg 


dS4  _  dm  (T  ^  ^  ^dEi  dQi 

wg-wg{Lhv+Ei)+mwg  +  wg~^{a9(P9~Pl)  +  Pi) 


ds4  _  .  dEr  dQi  ,  _  dvz 

dpi  m  dpi  dpi  d  d pi 


dS4  .  dEi  dFd  „dV,  dQi 
=  m—  +  +  Fd—  + 


r\  r\  1  '  l  r\  1  (in  1  r\ 

aui  oui  aui  dui  dui 


dS4  .  dEi  dQi 

m  =mwl  +  m~,i{(1~aa)(Pa~p,)~p,) 


(B.54) 

(B.55) 

(B.56) 

(B.57) 

(B.58) 

(B.59) 

(B.60) 


Rows  5,  6  and  7  are  the  opposite  sign  of  rows  2,  3  and  4  respectively. 
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APPENDIX  C: 

Time  Derivatives  and  Discretization  of 
Transversality  Condition  for  Single-  and 
Two-Phase  Control  Formulation 

For  the  single -phase  problem,  the  continous  form  of  the  functional  /  that  must  equal  zero  at 
the  optimal  final  time  is  given  in  Equation  2.98.  For  conveniance,  the  following  definition  is 
constructed,  suppressing  the  subscript  g  for  the  gas  phase. 

dP  dp 
dp  dt  + 
dP  dpu 
dpu  dt  ~*~ 
dP  dpE 
dpE  dt 
dP 
~dt 

The  discrete  form  is  shown  in  Equation  C.2. 

u(Xj,tN)2  p{Xj,tN)  -  p{xi,tN~1) 

2  At 

/  +N \  PU(Xi i  fN)  ~  PU(X^  ^_1)  , 
u{xi,r) - — - + 

pE(xh  tN)  -  pE(xj,  t^1) 

At 


fl 

dt 


bdP  dP 
dt  dt 


+  b(P  —  Q) 


d2P\ 

dt 2  J 


dx 


(C.3) 


The  continuous  derivative  of  the  term  is  given  in  Equation  C.4  and  the  discrete  form  in 
Equation  C.5. 
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du  dp  u2d2p  dudpu  d2pu  d2pE\  d2P 
dt  dt  +  2  dt 2  +  dt  dt  +  dt2  +  dt2  j  +  dt2 


=  (7  - 1) 


u(Xi,tN ) 


u(Xi,tN)  -  u(Xi,tN  *)  p(xi,tN)  -  p(xi,tN  *) 


”  '  At  At 

■u2(a7,  t^)  p(xj,  tN)  -  2p(xh  t^-1)  +  p(xi,  t^-2) 
2  At2 

-  u(xj,  t^1)  ptt(gjj,  t^)  -  pu(xj,  tN -1) 

At  At 


u(xi,tN) 


N.pu(xi,tN)-2pu(xi,tN  l)+pu(xi,tN  2) 


pE(xi,tN)  -  2pE(xi,tN  l)+pE(xi,tN  2) 
At2  + 

1  P(Xi,  tN )  -  2P(si,  t*"1)  +  Pfo,  tN~ 2) 
7  —  1  At2 


(C.5) 


df  ™(Xi,tN) 


z(Xi,tN)  -  z(Xi,tN  x)  P(Xi,tN)  -  P(Xi,tN  x) 

At  +  At 

^Pte^-Q^t*))  -±DPi 


DPi+ 


(C.6) 


For  the  two-phase  problem,  the  derivative  of  the  continuous  form  of  the  Transversality  Condi¬ 
tion  is  calculated  from  Equation  2.115  and  shown  here  in  Equation  C.7  and  it’s  discrete  sum¬ 
mation  form  in  Equation  C.8. 


df  r  (dpy  ^  d2p  r 

+{-p-®wdx 


(C.7) 
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$[_  ~  yM 

dt  ~  i=1 


' P{xhtN)  -  P^t^W 


\  (P(xi,tN)  -Q(xi)) 


At 

f-N 


+ 


P(xi,  tN)  -  2 P(Xi,  t"-1)  +  P(Xi,  tN -2) 
At2 


Ax 


(C.8) 
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